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ABSTRACT 


Modem computational and experimental tools for aerodynamics and propulsion applications have 
matured to a stage where they can provide substantial insight into engineering processes involving fluid 
flows, and can be fruitfully utilized to help improve the design of practical devices. In particular, rapid and 
continuous development in aerospace engineering demands that new design concepts be regularly proposed 
to meet goals for increased performance, robustness and safety while concurrently decreasing cost. To date, 
the majority of the effort in design optimization of fluid dynamics has relied on gradient-based search 
algorithms. Global optimization methods can utilize the information collected from various sources and by 
different tools. These methods offer multi-criterion optimization, handle the existence of multiple design 
points and trade-offs via insight into the entire design space, can easily perform tasks in parallel, and are 
often effective in filtering the noise intrinsic to numerical and experimental data. However, a successful 
application of the global optimization method needs to address issues related to data requirements with an 
increase in the number of design variables, and methods for predicting the model performance. In this 
article, we review recent progress made in establishing suitable global optimization techniques employing 
neural network and polynomial -based response surface methodologies. Issues addressed include techniques 
for construction of the response surface, design of experiment techniques for supplying information in an 
economical manner, optimization procedures and multi-level techniques, and assessment of relative 
performance between polynomials and neural networks. Examples drawn from wing aerodynamics, 
turbulent diffuser flows, gas-gas injectors, and supersonic turbines are employed to help demonstrate the 
issues involved in an engineering design context. Both the usefulness of the existing knowledge to aid 
current design practices and the need for future research are identified. 


Keywords: Global Optimization. Response Surface Methodology, Design of Experiments, Neural 
Networks. 
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1. INTRODUCTION AND SCOPE 

Modem computational and experimental fluid dynamics tools have matured to a stage where they 
can provide substantial insight into engineering processes involving fluid flows. This can help analyze the 
fluid physics as well as improve the design of practical devices. In particular, rapid and continuous 
development in the technology of fluid machinery demands that new design concepts be regularly proposed 
to meet goals for increased performance, robustness and safety while concurrently decreasing cost. 

Most aerospace system and component designs are conducted as open loop, feed-forward 
processes. For example, for rocket engines, currently, one design iteration for a given set of engine balance 
conditions takes up to several weeks with the blade geometry design sub-iteration phase taking several days 
each. The quest for an acceptable blade surface velocity distribution is accomplished with many ad hoc 
rules in what is really a manual trial-and-error process. A systematic approach capable of identifying 
optimum design and comparing possible trade-offs can significantly improve the productivity and shorten 
the design cycle. 

Objective and efficient evaluation of advanced designs can be facilitated by development and 
implementation of systematic optimization methods. To date, the majority of the effort in design 
optimization of fluid dynamics has relied on gradient-based search algorithms (Baysal and Eleshaky [I], 
Lambert et al.[27], Reuther et al. [54]). These methods work iteratively through a sequence of local sub- 
problems, which approximate objective and constraint lunctions tor a sub-region ot the design space, e.g., 
by linearization using computed sensitivities. Major challenges tor these optimization approaches are the 
robust and speedy computation of sensitivity coefficients (Elbanna and Carlson [10], Dadone et ai.[6]). 

Local optimization methods based on derivatives are commonly used in such engineering system 
design optimization problems (Sobieszczanski-Sobieski and Haftka [59]). On the other hand, global 
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optimization techniques also have been commonly used for engineering design optimization problems 
especially for multidisciplinary ones. In its current practice, the global design optimization method involves 
three primary steps (Figure 1 ): (a) Generation of individual data sets within the design space, (b) 
Interpolation among these data sets via some continuous functional representation; and (c) Optimization ot 
the objective function via a certain search strategy. Yet despite recent research advances, formal design 
optimization has yet to see practical use in real design scenarios. The reasons are tive-fold: 

(1) Engineering design, even within a single discipline, typically involves many parameters (and 
hence many degrees of freedom) rather than the handful demonstrated in most research papers. This 
renders unrestricted “brute force” search schemes too resource-intensive. 

(2) The objective functions are likely to be multi-modal or discontinuous over the broad design space, 
rendering gradient search methods insufficient by themselves. Furthermore, the usual practice to combine 
multiple goals into a single quantitative objective function is too restrictive. Qualitative goals are often 
required to correctly characterize a problem. (E.g., maximizing a turbine blade’s aerodynamic efficiency 
with a smooth, monotonic surface velocity distribution, while spreading heat load as uniformly as possible.) 
Furthermore, these goals may have arisen from diverse disciplines and are usually treated sequentially by 
different groups. 

(3) It is inadequate to think of the final product of a design process as a mere geometry. A “design” 
really encompasses a whole set of operating, manufacturing and project level decisions. 

(4) As the interaction between numerical simulation and physical test data becomes stronger, the 
future engineering knowledge base is likely to consist ot all sorts of heterogeneous data sources including 
test data, experimental data, past product experiences, semi-empirical modeling, and high fidelity 
simulations. Some data are anecdotal; others cover only small “patches” of the physical domain but are still 
useful for “reality checks”. A unified framework needs to be constructed for representation, capturing and 
mining of all these data types so the response functions can be continuously improved. 

With the above observations, global optimization methods are attractive because they have several 
advantages when compared to local gradient-based methods (Shyy et al. [58]): 

(1) They do not require calculation of the local sensitivity of each design variable, 

(2) They can utilize the information collected from various sources and by different tools. 


i 


(3) They offer multi-criterion optimization, 

(4) They can handle the existence of multiple design points and trade-offs, 

(5) They easily perform tasks in parallel, and 

(6) They can often effectively filter the noise intrinsic to numerical and experimental data. 

Among global approximation techniques, the Response Surface Methodology (RSM) has gained 
the most attention since it consists of a simple way of connecting codes from various disciplines 
(Sobieszczanski-Sobieski and Haftka [59]). The RSM is a collection of mathematical and statistical tools 
used in investigative experimentation by scientists and engineers (Bauer Jr. et al. [2]). The RSM approach 
replaces the objective and constraint functions by simple functions, often polynomials, which are fitted to 
the carefully selected points. Since RSM can utilize information collected from various sources and by 
different tools, it can also offer multi-criterion optimization, handle the existence of multiple design 
selections and related trade-offs, and address the noises intrinsic to numerical and experimental data. A 
main advantage of RSM is its robustness and intelligibility. Robustness and the smoothness of 
approximations in noisy environments are achieved by performing extra analyses, compared to the number 
of regression coefficients. This is a distinct advantage over derivative-based search algorithms, which may 
encounter difficulties in the presence of spurious local optima (Madsen et al. [32]). 

1.1 SCOPE 

In this article, we first review the basic concepts and methodologies, then assess the current statue, 
via case studies, of the global optimization techniques. Particular attention is paid to two different 
techniques to generate information to construct the response surface (RS) namely; Neural Network (NN) 
and polynomial-based Response Surface Methodology. Neural Networks are models that contain many 
simple linear and non-linear elements operating in parallel and connected in patterns (Greenman [14]). 
Polynomial -based RSM models the system with polynomials of assumed order and unknown coefficients. 
The solution for the set of coefficients that best fits the training data is a linear least square problem, 
making it trivial compared to the solution for NN. which involves a nonlinear training process. In this 
article, two neural network types, namely, Back-Propagation NN (BPNN) and Radial Basis NN (RBNN), 
are investigated. 
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The BPNN consists of multi-layer networks with differentiable activation function. The BPNN is 
the most employed NN type in the optimization literature (Carpenter and Barthelemy [4], Faller and 
Schereck [1 1], Fan et al. [11], Greenman [14], Greenman, and Roth [15] & [16], Huang et al. [20], Illi et al. 
[21], Kangas et al. [24], Lavretsky [28], Lawrence et al [29] & [30], Madavan et al.[31], Maghami and 
Sparks [34] & [35], Nikolaidis et al. [39], Norgaard et al. [40], Rai and Madavan [51], [52] & [53], Rosset 
al. [55], Sparks and Maghami [60], Stepniewski and Jorgenson [61], and Stepniewski et al. [62]). 

RBNN is a more recently developed multi layer network with a linear regression process that 
makes the mathematics simpler and computational costs lower (Orr [41], [42] & [43]). RBNN tends to have 
many more neurons than BPNN but can be configured faster for the same training data. The basic reason 
for this is that back-propagation neurons can have outputs over a large region of the input space, while 
radial-basis neurons respond to relatively small regions of the input space. Thus, larger input spaces require 
more radial-basis neurons for training. More detailed evaluation of RBNN and BPNN will be given in the 
following sections. 

Polynomial -based response surfaces and linear regression techniques have been originally 
developed to filter noise from experimental data. Sophisticated statistical tools are available for these 
purposes. One class of tools, design of experiments , is often used to select points for training that minimize 
the effect of noise on the fitted polynomial. A second set of tools, analysis of variance, is routinely used to 
identify polynomial coefficients that are not well characterized by the data, and are therefore overly 
sensitive to noise. Analysis of variance helps avoid overfitting of the data, which otherwise would result in 
the mapping of the noise. On the other hand, neural networks are much more flexible in functional torm, 
which means that they can be better suited to fit complex functions that are not easily approximated by 
polynomials. For example, when the physical system changes from one regime to another due to the 
presence of critical parameters, NN performs better than RSM. This advantage is particularly useful when 
there is very little numerical noise, and it is possible to obtain very accurate approximations to the 
underlying function (Papila et al.[49]). The relative strengths and weaknesses of NN and polynomial -based 
RSM are summarized in Table l. 

Table 2 summarizes the existing literature evaluating the relative performance of NN and 
polynomial-based RSM approximation. For example, Carpenter and Barthelemy [4] used NN and 
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polynomial -based approximations to develop RS for several test problems. It is demonstrated that two 
methods perform comparably based on the number of undetermined parameters. Rai and Madavan [51] 
investigated the feasibility of applying neural networks to the design of turbomachinery airfoils. Neural 
network approach is used for both function approximation and prediction. It is found that neural networks 
are quite efficient in both tasks. An aerodynamic design procedure that employs a strategy called 
parameter-based partitioning incorporating the advantages of both traditional RSM and NNs to create a 
composite response surface is described by Rai and Madavan [52] & [53]. It is shown that such method can 
handle design problems with higher dimensional problems than would be possible using NN alone. 
Nikolaidis et al. [39] used NNs and response surface polynomials to predict the performance characteristics 
of automotive joints using geometrical parameters. It is shown that both methods pertormed comparably. 
NN-based aerodynamic design procedure is applied to the redesign ot a transonic turbine stage to improve 
its unsteady aerodynamic performance by Madavan et al.[3l|. It is illustrated that using an optimization 
procedure combining the advantages of NN and polynomial-based RSM can be advantageous. Papila et al. 
[49] investigated the relative merits of polynomial -based RSM, RBNN and BPNN in handling dilferent 
data characteristics. It is demonstrated that using RBNN rather than BPNN has certain advantages as data 
size increases. Also, it is shown that RBNN gives more accurate results than polvnomi la-based RSM as the 
nature of the experimental data becomes complex. Shyy et al. [57] have integrated neural network 
techniques and polynomial -based RSM to obtain improved optimization tools. In Rai and Madavan [53], a 
composite NN & polynomial -based RS methodology is applied for a transonic turbine and it is 
demonstrated that a systematic application of such method can enhance the effectiveness of the overall 
optimization process. In the study by Vaidyanathan et al.[71 ], the application of NN and polynomial-based 
RSM in preliminary design of two rocket engine components, gas-gas injector and supersonic turbine, with 
modest amounts of data are discussed and it is demonstrated that NN and polynomial-based approximations 
can perform comparably for modest data sizes. 

In this article, we focus on the recent efforts in developing, improving, and optimizing appropriate 
techniques for design optimization of airfoils and rocket engine components capable of being used in 
applications like Reusable Launch Vehicles. Some of the physical components used as case studies are low 



Re aerodynamics, 2-D turbulent planar diffuser, the injector and the supersonic turbine for rocket 
propulsion. 

Specifically, the following issues are discussed: 

(1) The capability of the NN and polynomial-based RSM for handling data with variable sizes and noise, 

(2) The selection of NN configuration that is suitable for given design problems, 

(3) The effect of the design parameters on the performance of the NN, 

(4) The effect of distribution of the data over the design space in the construction of the global model, 

(5) The merit of employing a multi-level optimization strategy to perform the task adaptively and 
efficiently, 

(6) Possible trade-offs between capacity design objectives and their impact on design selections. 

2. REVIEW OF METHODOLOGIES 

In response-surface-based global optimization, there are several key technical elements: 

(1) Response surface with polynomials and statistical analysis 

(2) Neural networks with BPNN and RBNN 

(3) Design of experiments with face centered composite design (FCCD), orthogonal arrays (OA) and D- 
Optimal designs 

(4) Optimization procedure including the multilevel approach. 

In the following, we review these elements in sequence. 

2.1 RESPONSE SURFACE METHOD (RSM) 

The approach of RSM is to perform a series of experiments, based on numerical analyses or 
physical experiments, for a prescribed set of design points, and to construct a global approximation ot the 
measured quantity over the design space (Figure 1). The polynomial-based RSM, used in all the case 
studies referred to, constructs polynomials of assumed order and unknown coefficients based on regression 
analysis. The number of coefficients to be evaluated depends on the order of polynomial and the number of 
design parameters involved. For instance, a second-order polynomial of N design variables has 
(N+lXN+2y<2!) coefficients. A cubic model has (N+l )(N+2)(N+3)/(3!) coefficients. In this article, the 
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polynomial approximations are constructed by standard least square regression using JMP [23], a statistical 
analysis software that provides a variety of statistical analysis functions in an interactive manner. 

In the practical application of RSM, it is necessary to develop an approximate model for the true 
response surface. The approximate model is based on observed data from the process or system and is an 
empirical model. The second order (quadratic) response surface model for response variable v with k 
regressors can be written as 

y = a, + TP,*, + ypy +y.'L^r' c ' x j +£ < 1 > 

(=1 .=1 i=l /= 2 

The above equation can be written in matrix notation as follows 

y~XP+E ( 2 ) 

where y : ( nxl ) vector of observations, X: (nxn p ) matrix of the levels of the independent variables, p : ( n p xl ) 
vector of the regression coefficients, s: {nxl) vector of random error, n: the number of observations, and n p : 
the number of terms in the model. 

The purpose is to find the vector of least square estimators, b y that minimizes 
L = y e: = £ l E = (.y - X P) r {y -Xp) (3) 

; = l 

which yields to the least squares estimator of p 

b=(X ! XY l X r y < 4 ) 

The global fit and prediction accuracies of the response surfaces are assessed through statistical 
measures such as the t-statistic, or t-ratio, root mean square error (rms-error), variation (Myers and 
Montgomogery [37]). The t-statistic is determined by 

t=-^— ( 5 ) 

se{ b } ) 

where b t is the least squares estimators of / h regression coefficient and se(bj) is the standard error of bj and 
it is given by 

se(b l ) = <T ay JC~ ( 6 ) 
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where C n is the diagonal element of (X f X)~' corresponding to b } . Here < 7 a is the adjusted rms-error (or rms- 

error predictor) incurred while mapping the surface over the data set. The quality of the fit of the ditferent 
surfaces can be evaluated by comparing the adjusted rms-error value that is defined as: 



( 7 ) 


where e , is the error at i lh point of the training data. 

The accuracy of the models in representing the design space is gauged by comparing the values of 
the objective function at test design points, different from those used to generate the fit, with the empirical 
solution. The prediction rms-error, cr, for the test set is given by: 



In this equation e t is the error at the i th test point and m is the number of test points. 

The coefficient of multiple determination R 2 measures the proportion of the variation in the 
response around the mean that can be attributed to terms in the model rather than to random error and it is 
determined by 

R : = SS JL = [ _SS A _ ( 9 ) 

SS w SS VV 


SSf. : is the sum of squares of the residuals (=]T ( y, -y ,) 2 ) where y is the predicted value by the fitted 

r=I 

model. SSft is the sum of squares due to regression (=^(y, ~y)~ ) where y is the overall average ot v,. 

<=i 

SSyy is the total sum of squares about the mean given by 

SS w = SS h + SS, = yy y, - y f (10) 


where y is overall average of y v 

R 1 - is an R 2 value adjusted to account for the degrees of freedom in the model and is given by 


R : = [ _ SS f J(n-p) =l _ 


n - 1 


SS iv /(n- 1) \"~P 


(1 -/?“) 


(ID 
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Since R 2 increases as terms are added to the model, the overall assessment of the model may be better 
judged from R~ . 

The polynorrual -based response surface techniques are effective in representing the global 
characteristics of the design space. It can filter the noise associated with design data. Since, the solution for 
the set of coefficients that best fits the training data is a linear least square problem, it is trivial compared to 
the solution for the NN coefficients, which is often a non-linear least square problem. The linearity of the 
polynomial-based RSM also allows us to use statistical techniques known as design of experiments (DOE) 
to find efficient training sets. On the other hand, depending on the order of polynomial employed and the 
shape of the actual response surface, the RSM can introduce a substantial error in certain region of the 
design space. An optimization scheme requiring large amounts of data and a large evaluation time to 
generate meaningful results is hardly useful. 

2.2 NEURAL NETWORKS (NN) 

Neural networks are massively parallel computational systems comprised of simple nonlinear 
processing elements with adjustable interconnections. Neural networks simulate human functions such as 
learning from experience, generalizing from previous to new data, and abstracting essential characteristics 
from inputs containing irrelevant data (Greenman [14]). The processing ability of the network is stored in 
the inter-umt connection strengths or weights obtained by a process of adaptation to, or learning from, a set 
of training patterns. Training of a network requires repeated cycling through the data, each time adjusting 
the values of the weights and biases to improve performance. Each pass through the training data is called 
an epoch and the NN learns through the overall change in weights accumulating over many epochs. 
Training continues until the error target is met or until the maximum number of neurons is reached. In 
Figure 2, a neuron model with multiple inputs and bias is shown. 

Accordingly, the input is transmitted through a connection that multiplies it with the weight 
related to that connection. The bias is similar to a weight except that it has a constant input of l . The effect 
of the product weight and input and the bias are added at the summing junction to form the net input for the 
transfer (or activation) function. In Figure 3, a multi-layer network is shown. A layer of network includes 
the combination of weights, the multiplication and summing operations, the biases and the transfer 
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functions. In a layered neural network, neurons in every layer are associated with neurons in the previous 
layer in such a way that the outputs of an intermediate layer are the inputs to the following layer. The layer 
that produces the network output is called an output layer. All other layers are known as hidden layers. 

Even though research on neural network started in early 1940s, NN became quite popular around 
1980s with the introduction of multi-layered NN (Rumelhart et al. [56]) in a wide range of disciplines, 
including engineering. Over the last decade, NN approach has been used in the aerospace related industry. 
Illi et al. [21] examined the application of NN technology to an automated diagnostic and prognostic 
system for turbine engine maintenance. Preliminary results indicated that using NN to maintain diagnostics 
saves time and improves performance. Kangas et al. [24] used Back-Propagation NNs (BPNN) to monitor 
turbine engine performance and diagnose failures in real-time. The application of NN technology appears 
to hold great promise for enhancing the effectiveness of army maintenance practices. Huang et al.[20] 
developed and evaluated a multi-point inverse airfoil design method using NNs. It is shown that neural 
network predictions are acceptable for lift and moment coefficient predictions. Time dependent models that 
predict unsteady boundary layer development, separation, dynamic stall and reattachment are developed by 
Faller and Schereck [1 1] using NNs. It is demonstrated that NNs can be used to both predict and control 
unsteady aerodynamics effectively. Fan et al. [12] introduced a new approach for active laminar tlow 
control that incorporates BPNN into a smart wall interactive tlow control system. Convergence of the 
BPNN is investigated with respect to the complexity of the required function approximation, the size of the 
network in relation to the size of optimal solution and the degree of noise in the training data by Lawrence 
et al. [29]. The techniques and principles for the implementation of neural network simulators are also 
presented by Lawrence et al. [30]. Methods for ensuring the correctness of results avoiding duplication, 
automating common tasks, using assertions liberally, implementing reverse algorithms, employing multiple 
algorithms for the same task, and using extensive visualization are discussed. Efficiency concerns, 
including using appropriate granularity object-oriented programming, and pre-computing information 
whenever possible, are also studied. Norgaard et al. [40] used BPNN for more effective aerodynamic 
designs during wind tunnel testing. Four different NNs are trained to predict coefficients of lift, drag, 
moment of inertia, and lift drag ratio (Q, C D , Cm and UD) from angle of attack and flap settings. Hybrid 
neural network optimization method is successfully applied to produce fast and reliable predictions of 
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aerodynamic coefficients and to find optimal flap settings, and flap schedules. Ross et al. [55] applied 
BPNN to minimize the amount of data required to completely define the aerodynamic performance of a 
wind tunnel model. It is shown that the trained NN has a predictive accuracy equal to or better than the 
accuracy of the experimental measurements using only 50% of the data acquired during the wind tunnel 
test. BPNN is employed for rapid and efficient dynamics and control analysis of flexible systems by Sparks 
and Maghami [60]. It is demonstrated that NN can give very good approximations to nonlinear dynamic 
components, and by their judicious use in simulations, allow the analyst the potential to speed up the 
analysis process considerably once properly trained. The high-lift performance of a multi-element airfoil is 
optimized by using neural-net predictions by Greenman [14]. 

BPNN have been successfully integrated with a gradient-based optimizer to minimize the amount 
of data required to completely define the design space of a three-element airfoil. It is shown that using NN 
reduced the amount of computational time and resources needed in high-lift rigging optimization. 
Greenman and Roth [15] also applied BPNN for high-lift performance of a multi-element airfoil and it is 
demonstrated that the trained NN predicted the aerodynamic coefficients within an acceptable accuracy 
defined to be the experimental error. Stepniewski and Jorgenson [61] used a Singular Value Decomposition 
based node elimination technique and enhanced implementation of the Optimal Brain Surgeon algorithm to 
choose a proper NN architecture. It is demonstrated that combining both methods creates a powerful 
pruning scheme that can be used for tuning feed-forward connectionist models. Maghami and Sparks [34] 
& [35] also demonstrated that the methodology they developed based on statistical sampling theory 
guarantees that the trained networks provide a designer-specified degree of accuracy in mapping the 
functional relationship. BPNN is used to fill in a design space of computational data in order to optimize 
flap position for maximum lift for a multi element airfoil by Greenman and Roth [16]. A Genetic 
Algorithm (GA) and gradient-based optimizer are used together with NN and it is found that the 
demonstrated method has a higher fidelity and a reduction in CPU time when compared to an optimization 
procedure that excludes GA. Approximation abilities of BPNN is addressed by Lavretsky [28]. A novel 
matrix method for multi-input-multi-output NN is introduced and it is shown that by allowing inner layer 
connections as well as connections between any layers, ordered NN has superior interpolation ability when 
compared to conventional feed-forward NN. Stepniewski et al. [62 1 presented a new hybrid method that 
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combines a bootstrap technique and a collection of stochastic optimization method such as GA for 
designing a NN. The method minimizes generalization error. It is demonstrated that the solutions produced 
by this method improve the generalization ability on the average of five to six times when compared to 
pruned methods. 

All of the above listed references preferred to use BPNN among the other NN choices (Demuth 
and Beale[8], Kosko[25], and Jang et al.[22]). This is due to the fact that BPNN attempts to have small 
number of neurons when compared to the other NNs. However, since BPNN is usually slower because at 
each step the error is propagated back to all the weights in the system, other NNs could be more efficient 
than BPNN for specific problems. This article reviews the works focusing on Radial-Basis NN (RBNN) 
and BPNN models developed by using Matlab (Dernuth and Beale [8]). A comparative study for radial- 
basis and back-propagation approaches is also included. Brief summaries of the two approaches are given 
in the following sections. 

2.2.1 Back- propagation Networks (BPNN) 

Back-propagation networks are created by generalizing the Widrow-Hoff learning rule (Demuth 
and Beale [8] and Kosko[25]) to multiple-layer networks and nonlinear differentiable transfer functions. 
These networks are multi-layer networks with hidden layers of sigmoid transfer function and a linear output 
layer. The transfer function in the hidden layers should be differentiable and thus, either log-sigmoid or tan- 
sigmoid functions are commonly used. In this article, a single hidden layer with a tan-sigmoid transfer 
function, tansig (Figure 4), given as tanh(n), if n is the input is considered. The maximum and minimum 
outputs of the function are 1 and -1, respectively. 

The output of the function is given by 

a = tansig ( wx X +b) (12) 

where tansig is the transfer function, w is the weight vector, X is the input and b is the bias. For BPNN, the 
initial weights and biases are randomly generated and then the optimum weights and biases are evaluated 
through an iterative process. The weights and biases are updated by changing them in the direction of down 
slope with respect to the sum-squared error of the network, which is to be minimized. The sum-squared 
error is the sum of the squared error between the network prediction and the actual values of the output. In 
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BPNN (Figure 4a) the weights, w,, and biases, 6;, in the hidden tansig layer are not fixed as in the case of 
RBNN. Hence, the weights have a nonlinear relationship in the expression between the inputs and the 
outputs. This results in a nonlinear regression problem, which takes a longer time to solve than RBNN. 
Depending upon the initial weights and biases, the convergence to an optimal network design may or may 
not be achieved. Due to the randomness of the initial guesses, if one desires to mimic the process exactly 
for some purpose, it is impossible to re-train the network with the same accuracy or convergence unless the 
process is reinitiated exactly as before. The initial guess of the weights is a random process in Matlab. 
Hence to re-train the network the initial guess has to be recorded. 

The number of neurons in the hidden layer of a back-propagation network is a design parameter. It 
should be large enough to allow the network to map the functional relationship, but not too large to cause 
overfitting. As a rule of thumb to choose the number of neurons in the hidden layer, Greenman [14] used 
2s+l where s is the summation of total number of inputs and total number of outputs and Carpenter and 
Barthelemy [4| used m+1 where m is the number of nodes in the output layer. Once the number of neurons 
in the hidden layer is decided, the network design is reduced to adjusting the weighting coefficient matrices 
and the weighting bias vectors. These parameters for BPNN are usually adjusted using a gradient method 
such as the Levenberg-Marquardt technique (Greenman [14], Norgaard et al. [40], Ross et al. [55], Sparks 
and Mag h ami [60], Stepniewski et al. [62]). In Matlab, BPNN can be trained by using three different 
training functions, trainbp, trainbpx and trainlm. The first two are based on the steepest descent method. 
Simple back-propagation w'lth trainbp is usually slow since it requires small learning rates for stable 
learning. Trainbpx , applying momentum or adaptive learning rate, can be considerably faster method than 
trainbp but trainlm , applying Levenberg-Marquardt optimization, is the most efficient since it includes 
improvement techniques to increase the speed and reliability of simple back-propagation networks. The 
Levenberg-Marquardt update rule is 

AVV = (J r J + uI) l J r e 03) 

where SW is the change in woight, J is the Jacobian matrix of the derivatives of each error with respect to 
de 

each weight, i.e., — - , l is the identity matrix, fu is a scalar and e is the error vector. It the scalar p, is 
large, the above expression approximates the steepest descent, while if it is small then the method reduces 
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to the Gauss-Newton method. The Gauss-Newton method is faster and more accurate near an error 
minimum, so the aim is to shift towards the Gauss-Newton method as quickly as possible. Therefore, fi is 
decreased after each successful step and increased only when a step increases the error. The design 
parameters for trainlm are the number of neurons in the hidden layer, 5/, a user defined sum square error 
goal, and the maximum number of epochs. The training continues until either the error goal is reached, the 
minimum error gradient occurs, the maximum value of fu occurs, or the maximum number of epochs has 

been met. 

2.2.2 Radial-Basis Neural Networks (RBNN) 

Radial-basis neural networks are tw^o-Iayer networks with a hidden layer of radial-basis transfer 
function and a linear output layer. The main advantage of this approach is the ability of keeping the 
mathematics simple and computational costs low due to linear nature of RBNN (Orr [41]). Outline of 
supervised learning, main application area for RBNNs and the least squares method used together with 
supervised learning with linear models are explained in detail in Orr [41]. Optimum of the regularization 
parameter of RBNN is also searched in this paper. A computational method for re-estimating the 
regularization parameter of RBNN, based on generalized cross-validation, is explained by Orr [42]. RBNN 
is designed in such a way that it can adapt the width of the basis function, and it is found that it can predict 
better than a similar RBNN with the fixed width basis function. Orr [43] explains improvements made for 
to forward selection and ridge regression methods. A methodology that is a cross between regression trees 
and RBNN is described. The size of RBNN is also optimized based on regularization parameter in Orr [42]. 
The transfer function for radial-basis neuron is radbas , which is shown in Figure 5. Radbas , given 

as e~ n , where n is the input, has maximum and minimum outputs of 1 and 0 respectively. The output of 
the function is given by 

a = radbas[dist{ vv, X )x b) 0^) 

where radbas is the transfer function, dist is the vector distance between the network’s weight matrix, w, 
and the input vector, X and b is the bias. Radial-basis transfer function radbas calculates its output 

according to a - e"' 
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In a radial basis network (Figure 5a) each neuron in the radbas hidden layer is assigned weights, 
wj which are equal to the values of one of the training input design points. Therefore, each neuron acts as a 
detector for a different input. The bias for each neuron in that layer, b/ is set to 0.8326/5C, where sc is the 
spread constant, a value defined by the user. This defines the region of influence by each neuron. The 
training process is then reduced to the evaluation of the weights, w 2 , and biases, b 2 . in the output linear 
layer, which is a linear regression problem. If the input to a neuron is identical to the weight vector, the 
output of that neuron is 1, since the effective input to the transfer function is zero. When a value of 0.8326 
is passed through the transfer function the output is 0.5. For a vector distance equal to or less than 0.8326/6, 
the output is 0.5 or more. The spread constant defines the radius of the design space over which a neuron 
has a response of 0.5 or more. Small values of sc can result in poor response in a domain not closely 
located to neuron positions, that is, for inputs that are far from the training data as compared to the defined 
radius, the response from the neurons will be negligible. Large values will result in low sensitivity of 
neurons. Since the radius of sensitivity is large, neurons whose weights are different from the input values 
by a large amount will still have high output thereby resulting in a flat network. The best value of the 
spread constant for some test data can be found by comparing a for networks with different spread 
constants. 

In Matlab, radial-basis networks can be designed using two different design procedures, solverbe 
and solverb. Both procedures require a spread constant , sc, as a design parameter; i.e., the radius of the 
basis in the input space to which each neuron responds. Solverbe designs a network with zero error on the 
training vectors by creating as many radial-basis neurons as there are input vectors. Therefore, solverbe 
may result in a larger network than required and may fit the numerical noise. A more efficient design in 
terms of network size is obtained from solverb, which creates one neuron at a time to minimize the number 
of neurons required. At each epoch, neurons are added to the network until it satisfies a user specified error 
goal. The design parameters for solverb are the spread constant, error goal, and the maximum number of 
epochs whereas it is only the spread constant for solverbe. 

Radial-basts networks may require more neurons than a comparable BPNN. However, RBNN can 
be designed in a fraction of the time it takes to train the standard BPNN due to non-linear regression 
process of back-propagation networks. Therefore, RBNN are more efficient to train when there is a large 
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amount of training data available. In Papila et ah [49], an effort is made to compare the accuracy and 
computing requirements between the radial-basis and back-propagation approaches with different sizes ot 
training data. Vaidyanathan et al. [71] also investigated relative performances of RBNN and BPNN for gas- 
gas injector and supersonic turbine. As will be discussed in the following sections, it is illustrated that 
among all the NN configurations, RBNN designed with solverb seems to be more consistent in 
performance for different data sets and RBNN, even when designed efficiently with solverb, tend to have 
many more neurons than a comparable BPNN with tart-sigmoid or log-sigmoid neurons in the hidden layer. 
The basic reason for this is the fact that the sigmoid neurons can have outputs over a large region of the 
input space, while radial-basis neurons only respond to relatively small regions of the input space. 
Configuring a RBNN often takes less time than that required for a BPNN because the training process of 
RBNN is a linear in nature. 

2.3 DESIGN OF EXPERIMENTS (DOE) 

In RSM, selecting the representation of the design space is a critical step because it dictates the 
distribution of the information available for constructing the response surface. It is well established that the 
predictive capability of RSM is greatly influenced by the distribution of sampling points in design space 
(Unal et al. [69] & [70]). In order to select design points for training that minimizes the effect of noise on 
the fitted polynomial, design of experiment ( DOE) techniques can be applied. There are different types ot 
design of experiments techniques in the literature as reported by Haftka et al. [19]. For example, Unal et al. 
[70] discussed the D-optimal design for the representation of the design space for a wing-body 
configuration of a launch vehicle. They showed that D-optimal design provides an efficient approach for 
approximating model building and multidisciplinary optimization. Papila and Haftka [48] also applied face 
centered composite design ( FCCD ) to select the experiment points in the design space when approximating 
wing structural weight. Unal et al. [68] & [69] studied response surface modeling using orthogonal arrays 
(OA) in computer experiments for reusable launch vehicles and illustrated that using this technique can 
minimize design, development, test and evaluation cost. Unal and Dean [67] studied the robust design 
method based on the Taguchi method (Unal and Dean [66] and Dean [7]) to determine the optimum 
configuration of design parameters for performance, quality and cost. They demonstrated that using such a 
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robust design method for selection of design points is a systematic and efficient approach for determining 
the optimum configuration. Brief summaries of FCCD , OA , and D-Optimal designs are given below. 

2.3.1 Face Centered Cubic Design (FCCD) 

Face centered cubic design (FCCD) creates a design space composed of eight corners of a cube, 
centers of the six faces and the center of the cube. Figure 6 shows FCCD selections for three design 
variables. The FCCD yields (2 ' v +2jV+ 7) points, where /V is the number of design variables. It is more 
effective when the number of design variables is modest, say, not more than 5 or 6. The FCCD is used tor 
fitting second-order response surface. 

2.3.2 Orthogonal Arrays (OA) 

An orthogonal array (OA) is a fractional factorial matrix that assures a balanced comparison of 
levels of any factor or interaction of factors. Consider A, a matrix with elements of a\ where; shows the 
row (j = 1,2... n r ) and i shows the column (/ = 1,2 ...nc) that a[ belongs to, supposing that each a] e Q = 
{0, 1 . ..q -[ }. A is called an orthogonal array of strength t <n t if in each n r -row-by-r-column sub-matrix of 
all q l possible distinct rows occur A times (Owen [44]). Such an array is denoted by OA(n r ,n r ,q.t) by Owen 
[44]. 

Since the points are not necessarily at the vertices, the OA can be more robust than the FCCD in 
interior design space and are less likely to fail the analysis tool. Based on the design of experiment theory, 
orthogonal arrays can significantly reduce the number of experimental configurations. 

2.3.3 D-Optimal Design 

A D-Optimal design minimizes the generalized variance of the estimates, which is equivalent to 
maximizing the determinant of the moment matrix, V/ (Myers and Montgomery [37]). 



where n is the number of observations and n p is the number of terms in the model. 

The D-Opiimal Design approach makes use of the knowledge of the properties of polynomial 
model in selecting the design points. This criterion tends to emphasize the parameters with the highest 
sensitivity (Haftka et al. [19]). 

2.4 OPTIMIZATION PROCESS 

2.4.1 Search Procedure 

The entire optimization process can be divided into two parts: (1) RS/NN phase for establishing an 
approximation, and (2) Optimizer phase. 

In the first phase, polynomials or NN models are generated with the available training data set. In 
the second phase the optimizer uses the RS/NN during the search for the optimum until the final converged 
solution is obtained The initial set of design variables is randomly selected from within the design space. 
The flowchart of the process is shown in Figure 7. 

The optimization problem at hand can be formulated as min[f{x )\ subject to lb < x < ub , where lb 
is the lower boundary vector and ub is the upper boundary vector of the design variables vector x If the 
goal is to maximize the objective function then fix) can be written as -g(x), where g(x) is the objective 
function. Additional linear or nonlinear constraints can be incorporated if required. The optimization 
toolbox in Matlab is used here employs a sequential quadratic-programming algorithm. 

2.4.2 Objective Function 

When attempting to optimize two or more different objective functions, conflicts between them 
arise because of the different relationships they have with the independent parameters. An equation 
expressing the relationship between opposing effects of performance and weight can be employed as a 
criterion to guide the optimization task. Both NN and polynomial-based RS techniques can handle such 
multi-criteria optimization tasks in a straightforward manner by building a composite response surface 
from individual response surfaces. Such a task would have been impossible without response surface. This 
composite response surface is referred to as the desirability function. The maximization of the composite 
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function effectively provides a compromise between the individual functions. An average of some form is 
normally used to represent the composite function A geometric mean is a solution, which gives a 
composite function of the form: 

e i y 

£> = 104 I (16) 

where D is the composite objective function, d t s are normalized values of the objective functions and / is 
the number of objective functions. Each of the d l are weighted depending upon the importance of the 
specific objective function. Figure 8 shows a typical trend for a desirability function with respect to the 
weighting factors. 

Another way of constructing a composite function is to use a weighted sum of the objective 
functions. The composite function can then be expressed as: 

D = Y^aJ, (17) 

i=l 

where D is the composite objective function and f 's are the non-normalized objective functions. The a,' s 
are dimensional parameters that control the importance of each objective function. 

3. DESCRIPTION OF THE CASE STUDIES 
3.1 GAS-GAS INJECTOR ELEMENT FOR ROCKET PROPULSION 

Development of an optimization scheme for injector design called method i ( Meth odology for 
Optimizing the Design of_Jnjectors) has been reported by Tucker et al. [64] & [65] Method i is used to 
generate appropriate injector design data and then guide the designer toward an optimum design subject to 
his specified constraints. As reported, method i uses the polynomial based RSM to facilitate the 
optimization. The RSM approach is to conduct a series of well-chosen experiments (empirical, numerical, 
physical or some combination of the three) and use the resulting information to construct a global 
approximation (response surface) of the measured quantity (response) over the design space. A standard 
constrained optimization algorithm is then used to interrogate the response surface for an optimum design. 
Neural network was also used in the design of shear co-axial injector element by Shyy et al. [57], and 
Tucker et al [64) & [65] along with the polynomial-based RSM.. Three different injector types are 
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considered, namely, shear co-axial injector element, impinging injector element, and swirl co-axial injector 
element. 


3.1.1 Shear Co-axial Injector Element 

The initial demonstration of method i by Tucker et al. [64] focused on a simple optimization of a 
shear co-axial injector element (Figure 9) with gaseous oxygen (GOj and gaseous hydrogen (GHj 

propellants. The goal is to maximize the energy release efficiency, ERE while minimizing the chamber wall 
heat flux, Q. This is achieved by maximizing a composite objective function given by: 

D = {d, :l>F d Q f 08) 

where the normalized functions are defined as in Eqns (19) and (20). In the case where a response should 
be maximized, such as ERE , the normalized function takes the form: 

d = f for A < ERE < B (19) 

” ( B-A ) 

where B is the target value and A is the lowest acceptable value. d FRi: is set to / for any ERE > B and d FHK - 
0 for ERE < A. The choice of s is made based on the subjective importance of this objective in the 
composite desirability function In the case where a response is to be minimized, such as Q , the normalized 
function takes on the form: 

(or C<Q<E (20) 

where C is the target value and E is the highest acceptable value. d,j is set to 1 for any Q < C and d L , = 0 for 
Q > E. A, B, C, and E are chosen according to the designer’s priorities or, as in the present article, simply 

as the boundary values of the domain of ERE and Q The value of t is again chosen to retlect the 

importance of the objectives in the design. In the study carried out, A and B are equal to 95.0 and 99.9, 
respectively. Values of C and E are equal to 0.48 and l.l, respectively. 

The design data was generated using an empirical design methodology developed by Calhoon et 
al. [3]. These researchers conducted a large number of cold-flow and hot-fire tests over a range of 
propellant mixture ratios, propellant velocity ratios and chamber pressure for shear co-axial, swirl co-axial. 
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impinging, and premixed elements. The data were correlated directly with injector/chamber design 
parameters, which are recognized from both theoretical and empirical standpoints as the controlling 
variables. For the shear co-axial element, performance, as measured by energy release efficiency, ERE , is 
obtained using correlations taking into account combustor length, L comb (length from injector to throat) and 
the propellant velocity ratio, V/V . The nominal chamber wall heat flux at a point just downstream of the 

/ o 

injector, Q nom , is calculated using a modified Bartz equation and is correlated with propellant mixture ratio, 
O/F, and propellant velocity ratio, V/V to yield the actual chamber wall heat flux, Q. The objective in the 

/ o 

initial demonstration of method i was to maximize injector performance while minimizing chamber wall 
heat flux (lower heat fluxes reduce cooling requirements and increase chamber life) and chamber length 
(shorter chambers lower engine weight). The data used to generate the polynomials and train the network is 
given in Table A3-A5. The quality of the response surface and neural networks are decided using 20 
additional design points different from those used to generate the models (Table A6). 


3.1.2 Impinging Injector Element 

The empirical design methodology of Calhoon et al . [3 ] uses the oxidizer pressure drop, AP 0 , fuel 
pressure drop, J P f , combustor length, and the impingement half-angle, a as independent variables. 

For this injector design, the pressure drop range is set to 10-20fr of the chamber pressure due to stability 
considerations. The combustor length, defined as the distance from the injector to the end of the barrel 
portion of the chamber, ranges from 2-8 inches. The impingement half angle is allowed to vary from 15- 
50°. Dependent variables include ERE (a measure of element performance), wall heat flux, <2^ injector heat 
flux, Q inp relative combustor weight, W reh and relative injector cost, C re i- 

The conditions selected for this example are: 

F = 1000 p si 

MR = 6 

( 21 ) 

m (W - 0.25 lb m / sec 
m GH =0. 0421b m sec 

The gaseous propellants are injected at a temperature of 540R. As noted above, the empirical 
design methodology used to characterize the ERE and Q w was developed by Calhoon et al. [3]. This 


methodology uses a quantity called the normalized injection momentum ratio to correlate the mixing at the 
different design points for the triplet element. They define this quantity as 


=' 


2.3mu 
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The maximum mixing, and thus maximum ERE , occurs at an MR nl of 2.0. Since the propellant 
mass flow rates are fixed, only the propellant velocities and the impingement half-angle influence the 
normalized injection momentum ratio. The velocities are proportional to the square root of the respective 
pressure drops across the injector, AP t , and A P } . For the flow conditions and variable ranges considered in 
this problem, MR ni ranges from 3.2 to 17.8. Accordingly, lowering AP 0 , raising AP f , increasing a, or some 
combination of these actions will increase ERE. The wall heat flux is correlated with the propellant 
momentum ratio as defined by 
mu 


MRt 


(23) 


For the F-O-F triplet element, i.e. the impingement injector element, the maximum wall heat flux 
occurs at a momentum ratio of approximately 0.4. High heat flux is the result of over-penetration of the 
fuel jet, which produces a high O/F in the wall region. For the flow conditions and variable ranges 
considered in this effort, MR ranges from 1.06 to 2. 1 1. Hence, increasing the value of this ratio by either 
increasing AP (} or decreasing AP, lowers the wall heat flux. 

The heat flux seen by the injector face. Q n ., is qualitatively modeled by the impingement height, 
H impmKe . The notion being that, as the impingement height decreases, the combustion occurs closer to the 
injector face, causing a proportional increase in Q inj . Thus, for the purposes of this exercise, Q lrxJ is modeled 
as the reciprocal of the H impinKe . Impingement height is a function of a and AP f . Reference to Figure 10 
shows that as a; is increased, is shortened. The dependence of H imptnKt on the fuel orifice diameter, 

df, and thus, AP, % results from making the freestream length of the fuel jet, L 1S , a function of d } . For each 
AP h L fs was set to six times d f for an impingement half-angle of 30°. So, as d, increases (corresponding to 
decreasing AP f ), L (s increases, as does 

The models for W rel and C re( are simple but represent the correct trends. W rr , is a function only of 
L i0mh , the combustor length from injector face to the end of the chamber barrel section. The dimensions of 


the rest of the thrust chamber assembly are assumed to be fixed. So, as L comb increases, W rel increases 
accordingly. The model for C re i is based on the notion that smaller orifices are more expensive to machine. 
Therefore, C rei is a function of both propellant pressure drops. As the AP increases, the propellant velocity 
through the injector increases and the orifice area decreases. So, as either, or both, AP 0 and AP t increase, 
Crei increases. 

The system variables given above and independent variables (constrained to the previously noted 
ranges) are used to generate the design data for element optimization studies. Since propellant momentum 
ratio is an important variable in the empirical design methodology, a matrix of momentum ratios was 
developed over the 100-200psi propellant pressure drop range. The matrix of 49 combinations of fuel and 
oxidizer pressure drops is shown in Table A7 where momentum ratios range from 1.06 to 2.11. Nine 
pressure drop combinations, eight around the border and one in the middle, were selected for use in 
populating the design database. These nine points are highlighted in Table A7 in bold type. 

Detailed design results for the case with both AP a and AP/ at 200psi are shown in Table A8. 
Similar data was generated for the other eight pressure drop combinations. There are 20 combinations ot 
Lcomb and a for each AP combination, making a total of 180 design points selected. Seventeen of these were 
outside the database embodied by the empirical design methodology, resulting in 163 design points actually 
being evaluated. The data trends are as expected. ERE , for a given AP combination, increases with 
increasing L comb and a. The increased L a)mh provides more residence time for the propellants to mix and 
bum. Increasing a increases the radial component of the injected fuel, thus providing better mixing. The 
wall heat flux is constant for a given AP combination. Impingement height increases with increasing a 
Relative combustor cost increases with increasing Lc 0mb and the relative injector cost is constant for a given 
AP combination. 

3.1.3 Swirl Co-axial Injector Element 

The chamber pressure, mixture ratio, and propellant How rates selected for this example are: 

P - 1000 psi 

MR = 6 

m :rk =025lb m ( sec 

m . M - 00421b ^ sec 
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The gaseous propellants are injected at a temperature of 540 R. Reference to Figure 1 1 shows that 
the G0 2t flowing in the center post of the element, exits the element with both radial and axial velocity 
components. This effect is achieved by introducing the G0 2 tangentially into the center post through small 
slots. When the G0 2 , under hydrostatic head, is forced through the tangential slots, part ol the pressure 
head is converted into a velocity head, causing a rotational velocity in the element. For a specified AP 0 and 
swirl angle, <9, the number and size of tangential slots, the discharge coefficient, the G0 2 center post 
diameter, d 0 , and the radial and axial G0 2 velocity components, V or and V oa are calculated. These quantities 
are then used to determine the dependent variables for each design condition. 

The element ERE , calculated according to the empirical design methodology of Calhoon et al.[3], 
is a function of all four independent variables noted above. A cold flow mixing efficiency, E rtl , w , for <9=90°, 
is correlated by: 


^.90 


100-5 In 



(25) 


The cold flow mixing length, L co{dy is correlated from a known chamber length, L cnmb . The G0 2 
post diameter, d 0 , is a function of AP„ and © Smaller values of d correspond to large values of AP t> and 
smaller swirl angles. The empirical swirl factor, is a function of the normalized differential injection 
velocity, (V r V„yV 0 . K s increases with increasing normalized differential injection velocity for the range of 
propellant velocities considered in this effort. For fixed propellant mass flow rates, the velocities V 0 and V } 
are functions of their pressure drops across the injector, AP () and AP fy respectively. For a given AP { „ V 0 also 
depends on the swirl angle. Lower V 0 's are a product of higher swirl angles. Cold flow mixing is thereby 
enhanced with higher values of V a (i.e. AP I} ) and L amh . Lower values of V } (i.e. AP ) and (9 also tend to 
enhance cold flow mixing. 


A fractional factor,/, is applied to E m yo to account for the lower levels of cold flow mixing found 
with swirl angles less than 90°. The resultant measure of cold flow mixing, E m , & is a product of E m 9il and 
/. This factor, for a given design, is a function of the normalized differential injection velocity and the ratio 
of radial to axial G0 2 velocity, V r f ,/V f „ r Increasing values of both quantities increase /, with a value of/ =1 
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being found at V ot /V oa =1 (<9=90°) for all values of (V/ - V o y V 0 . Larger values of f s increase cold flow 
mixing. These values are found at low AP 0 and high AP f and <9. There is no dependency of f s on chamber 
length. These trends are opposite those noted above. Finally, ERE is proportional to E m , & 


The wall heat flux Q w , is correlated with the propellant momentum ratio as defined by: 


MR = 


(26) 


The wall heat flux curve from the Calhoon et al. [3] methodology is fairly flat, varying only about 
10% from high to low for the range of pressure drops considered in this effort. Q w decreases with 
increasing V 0 (high AP a and low G) and decreasing V f (low APf). That Q w would decrease with increasing 
V 0 is counter to intuition. It seems that high values of V 0 , for any <9, would result in higher mixture ratios in 
the wall region, as is the case for liquid O 2 . Calhoon et al. [3] do not discuss this effect. 

The heat flux seen by the injector, Q inp is actually modeled by the distance from the injector at 
which the propellant streams intersect. This axial distance is measured at the radial position corresponding 
to the center of the co-axial fuel annulus, or gap. It is here that the streams begin to mix and bum. This 
measure is qualitative, but captures the trend that higher injector heat fluxes occur the nearer the injector 
that the combustion begins. The axial distance is affected directly by the swirl angle, and indirectly by the 
propellant pressure drops. Q inJ decreases with decreasing swirl angle, increasing GO 2 pressure drop and 
decreasing GH 2 pressure drop. Swirl angle has the largest effect, while AP ( , is the least significant factor. 

The relative combustor weight, W reh is simply a function of the combustor length, L comht the 
distance from the injector to the end of the barrel portion of the chamber. The longer the combustor, the 
more it weighs. 

The relative injector cost, C r ,/, is a function of the fuel gap width and the width of the tangential 
slots used to induce the swirl in the G0 2 center post. Larger values of both variables result in lower 
machining costs, and thus lead to lower injector cost. The fuel gap width increases with increasing AP l)t and 
decreasing values of AP ; and ©. Swirl slot width increases with lower values of AP (J and G. Overall, C re i 
decreases with increasing AP (> and decreasing APf and G. Fuel pressure drop and swirl angle are the most 
significant factors. 
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A matrix of propellant pressure drop combinations was developed and nine combinations were 
selected for use in populating the design database. There are 20 combinations of L loml> and G tor each A P 
combination, making a total of 180 design points selected. 

In the work by Tucker et al. [641 & [65], method i uses the Response Surface Method (RSM) to 
find optimal values of ERE , Q fnp W rt!l and C re \ for acceptable values of AP () , AP, % L romb and <9. The 

approach of RSM is to perform a series of experiments, or numerical analyses, for a prescribed set of 
design points, and to construct a response surface of the measured quantity over the design space. In the 
present context, the five responses of interest are ERE , Q w , Q inp W ret and C reh The design space for each 
element consists of the set of relevant design variables AP 0 , AP f , L comb and & The response surfaces are fit 
by standard least-squares regression with a quadratic polynomial using JMP [23]. A backward elimination 
procedure based on t-statistics is used to discard terms and improve the prediction accuracy. Five full 
quadratic response surfaces are constructed by using JMP. 

In the current case, it is desirable to maximize ERE and while simultaneously minimizing Q w , Q inp 
W rei and C rel . 

3.2 SUPERSONIC TURBINE FOR REUSABLE LAUNCH VEHICLES 

Supersonic turbines that drive fuel or oxidizer turbopumps in rocket engines are of great interest to 
the next generation space propulsion industry, including the Reusable Launch Vehicles (RL\ ). They are 
complex, high-speed devices that produce shaft power by ducting the flow ot hot gasses over specially 
shaped blades on a wheel. For rocket engine applications, maximizing the vehicle payload for a given 
turbine operating condition is the ultimate goal. The flow path should be designed in such a way that it 
wastes less energy so that turbine temperatures or the mass flow rate can be reduced, or the turbine can be 
made smaller, increasing the efficiency (or specific impulse) of the rocket engine. Any gain in turbine 
efficiency will be reflected in reduced propellant consumption, resulting in an increase in the payload. 
However, higher turbine performance usually entails multistage designs, which are heavier. The design of a 
supersonic turbine often involves a considerable number of design variables with structural and 
aerodynamic constraints. With the number of design parameters involved, the overall procedure of design 
optimization of supersonic turbines becomes a challenging task. 
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Papila et al. [501 have conducted a global optimization investigation to perform the preliminary 
design of the supersonic turbines, including the selection of the number of stages and design variables. 
From 1- to 2- to 3-stage turbines, the number of design variables increases substantially. In shape design, 
from vane to blade, from stage to stage, and from 2-D to 3-D, not only does the number of design variables 
increase, but also the interactions among design variables become more complicated. Papila et al. [50] 
intended to investigate the individual, as well as collective effects of design variables by varying the design 
scope systematically. Vaidyanathan et al. [71] have used the data of the 1 -stage turbine to conduct a 
comparative study between RSM and NN. 

For preliminary design stage, single-, two- and three-stage turbines are considered. The design 
variables can be separated into two categories, one related to geometry and the other to performance. They 
are summarized as fol lows: 

(1) Geometric inputs: The geometric inputs are needed to layout the turbine meridional geometry, e.g., 
mean diameter, last rotor annulus area, blade height ratio between the I st vane and the last rotor blade 
(linear distribution of blade heights is assumed between the 1 st vane and the last rotor blade), vane and 
blade axial chords. 

(2) Performance inputs: The performance inputs are needed to calculate the turbine efficiency, e.g., speed 
(RPM), number of stages, blade row reaction, and work split (if more than l stage is investigated). 

For single-stage turbine, 6 design parameters (Table 3) are selected. These are (1) the mean 
diameter, (2) speed ( RPM ), (3) exit blade annulus area, (4) vane axial chord, (5) blade axial chord, (6) stage 
reaction. 

For 2-stage turbine, there are 1 1 design parameters (Table 3), namely, (1) the mean diameter, (2) 
RPM, (3) exit blade annulus area, (4) 1st blade height (% of exit blade), (5) I st vane axial chord, (6) 1 st 
blade axial chord, (7) 2 nd vane axial chord, (8) 2 lld blade axial chord, (9) 1 st stage reaction, (10) 2 nd stage 
reaction, and (11) 1 st stage work fraction. Note that 2 nd stage work traction is not a design parameter since it 
can be calculated by using 1 st stage work fractions, i.e., w j2 =l-Wfi. 

There are 15 (Table 3) design parameters for 3-stage turbine. These are (1) mean diameter, (2) 
speed (RPM), (3) exit blade annulus area, (4) 1 SI blade height (% of exit blade), (5) 1 st vane axial chord, (6) 
I st blade axial chord, (7) 2 nd vane axial chord, (8) 2 nd blade axial chord. (9) 3 rd vane axial chord, (10) 3 rd 
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blade axial chord, (11) I st stage reaction, (12) 2 nd stage reaction, (13) 3 rd stage reaction, (14) 1 st stage work 
fraction, (15) 2 nd stage work fraction. Note that 3 rd stage work fraction is not a design parameter since it can 
be calculated by using 1 st and 2 nd stage work fractions, i.e., vv^=/-(wyy+vy, 2 ,). 

The composite objective function chosen by Papila et al. [50] for design optimization corresponds 
to the payload increment, Apay, versus turbopump efficiency and weight. The relation between Apay and 
these two parameters can be developed as follows based on mission prolile studies, engine balance 
perturbation and some detailed turbopump layout and stress information gained from other proprietary 
programs. 

Apay =CiX(Ti-r| b ) xlOO-(W-W b ) (27) 

where T| b is the baseline efficiency and W b is the baseline weight. The constant c ; indicates that for every 
point in efficiency gained, the amount of payload capacity of the RLV is increased cj per turbopump 
Therefore, Apay function represents the amount of increase in payload capacity. The results of both payload 
increment based and composite desirability function based optimization are illustrated for 1, 2, and 3-stage 
designs. The results of both payload increment based- and composite desirability function-based 
optimization are illustrated for 1, 2, and 3-stage designs in the following chapters. 

Two structural constraints are considered by Papila et al. [50]. In axial turbines the product ot the 
blade exit annulus area and the RPM square, i.e., AN 2 is an indication of the blade centrifugal stress, which 
should bind the speed of the turbine. In addition, the disk stresses are also a restriction. In turbomachinery 
industry, the maximum stress value due to disk burst is often represented by a pitchline velocity limit, i.e., 
V pilch . The pitchline velocity can be calculated by multiplying RPM and the mean radius. 

3.3 TURBULENT PLANAR DIFFUSER 

The goal was to accomplish maximum pressure recovery by optimizing the wall contours. The 
flow' is incompressible and fully turbulent with a Reynolds number of 10 5 , based on the inlet throat half- 
width, D. The overall geometry is defined by the ratio of inlet and outlet areas, and the diffuser length to 
height ratio. In this study the length to height ratio is fixed at 3.0, and the area ratio at 2.0. The shape of the 
diffuser wall is designed for optimum performance, with five design variables represented by B-splines. 
The CFD model is based on the full Reynolds-averaged Navier-Stokes equations, w ith the k-e two-equation 
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turbulence model in closure form. At the inlet of the flow domain, a uniform flow distribution is specified. 
Detailed discussion of this study can be found in Madsen et al. [33]. 


3.3.1 Objective 

The dimensionless pressure recovery coefficient C p is introduced as the objective function to be 
maximized 

F = c = .. 4 g _ (28) 

' 1 , 

Here Ap is the static pressure difference between channel cross sections up- and downstream of 
the diffuser respectively, p is the fluid density, and u !nt is the inlet mean velocity. Inlet and outlet static 
pressures are averaged, even though the pressure distribution is nearly uniform due to well-developed flow 
at the considered cross sections. The CFD model uses a symmetry condition along the channel center axis, 
and has a computational mesh consisting of 120x50 cells including a long outlet section to establish a fully 
developed exit profile. The overall geometry of the two-dimensional planar diffuser, see Figure 12, is 
defined by the ratio of inlet and outlet areas, AR, and the diffuser length/height ratio, L=D, where L is the 
axial length of the diffuser. In this study the D=D-ratio fixed at 3 0, and the area ratio AR at 2.0. Expressed 
in terms of the inlet half-width D, the horizontal position of the inlet is ID, while the horizontal position of 
the outlet is 10D. The shape of the diffuser wall is designed for optimum performance, and to this end two 
separate cases of wall parameterizations are tried: (I) a two design variable case, where a polynomial 
describes wall shapes, and (2) a five design variable case that uses B-splines. Even though two different 
curve descriptions are used in the two cases, the most noteworthy difference seen from the point of view of 
the RSM lies in the problem size. 


3.3.2 Geometric Representation 

For shape parameterization in more variables, B-splines were preferred to natural splines 
(piecewise polynomials), although the latter technique is closer to the polynomial representation. B-splines 
excel in the predictable way that control points influence curve shape, and in the local control , which 
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prevents small changes in a control point position from propagating over the entire curve. Combined with 
low computational cost, these advantages have contributed to B-spline curves becoming a standard 
geometric modeling technique in computer-aided design. 

A B-spline is given in a parametric form as p(u): 

^ (,<)= PT!l = a9) 

_v(m)J is 

A set of blending functions N t:k combines the influence of n+1 control points P t , over the range of 
the parametric variable u. The blending functions N l;k are recursively determined polynomials with degree 
k-U where the parameter k dictates the order of continuity of the curve, and thus how many control points 
influence a curve segment. In this work k is 8, which corresponds to C 6 -continuity. The number of control 
points is 8 as well-two endpoints, five design variables and one point used for prescribing the inlet slope. 

B-splines have an approximating nature, in that they do not necessarily pass through control 
points, except for fixed curve endpoints. The slope at a curve endpoint is tangential to a straight line 
connecting the endpoint and the first control point, and may be prescribed by placing an additional fixed 
control point near the endpoint. 

Experimental and numerical evidence indicates that maximum pressure recovery in diffusers 
occurs at the border of appreciable flow separation. For this reason, strongly separated diffuser flows 
should be avoided, which makes it reasonable to restrict the design space to monotonic wall shapes. While 
the approximation accuracy does of course benefit from the reasonable design space approach, it is equally 
important in the present example that monotonicity constraints eliminate convergence problems associated 
with CFD-analysis of odd, non-monotonic designs. 

The parametric form in which B-splines are defined makes it non-trivial to derive monotonicity 
constraints analytically, so instead a constraint approximation G was set up in the form of a response 
surface for the minimum wall slope G. Then, observing the inequality constraint G >0 implies a positive 
wall slope and thus monotonicity throughout. Since B-splines are inexpensive to generate, 9 5 (59049) B- 
splines were computed (this took just a few seconds) and used for fitting a quadratic response surface. The 
approximation to the monotonicity constraint precludes some designs that satisfy the exact monotonicity 
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requirement. However, the effect of these inaccuracies on the solution of the optimum design problem is 
negligible. 

The regression analysis, to find 21 polynomial coefficients in five dimensions, is based on a 35- 
point D-optimal design. The surplus of analyses is generally required for reducing the sensitivity to 
numerical noise and to errors due to the simplified representation as a quadratic polynomial. Again, a pool 
of candidate points was created, this time using nine levels for each variable, (values ranging from 0.0 to 
1.0), and then checking the monotonicity of the B-splines for each of the 9 5 = 59049 designs. Observe that 
limiting the y-coordinate of the control points to a variation in the range [0:0; 1:0] is a somewhat artificial 
requirement, as monotonic shapes exist with coordinates slightly outside this range. A total of 20864 points 
are monotonic in wall shape. This relatively large percentage of acceptable cases reflects the smoother 
nature of approximating curves. Had a non-segmented polynomial curve representation been used, the 
condition of monotonicity in the control points would alone have reduced the number ot teasible design 
points to less than 1% of those inside a five dimensional box. As in the two-design-variable case, the subset 
of D-optimal points was found using the JMP- soft ware. 

3.4 LOW REYNOLDS NUMBER WING MODEL 
3.4.1 Training Data 

The aerodynamic model, a rectangular wing with a NACA 5405 airfoil cross-section (Figure 13) is 
designed for low Reynolds Number (Re=I(f-l(f) flows. Since airfoil performance decreases at low 
Reynolds Number flights, attempts to shrink the overall aircraft size while trying to keep sufficient lifting 
areas, result in low aspect ratio wing plan forms. As aspect ratio decreases, the percentage of the wing area 
affected by the tip vortices increases, creating a 3-D flow field over most of the field. Therefore, the 
analysis of such flows should consider the effects on performance and the effects of both the airfoil 
geometry (such as maximum camber) and the wing geometry (such as aspect ratio). In this study, the 
aerodynamic analysis is based on a 3-D potential flow solver, PMARC , and a 2-D coupled inviscid-viscous 
flow solver, XFOIL . The lift coefficients, Q, and drag coefficients, C 0 , for various maximum camber, y c , 
aspect ratios, AR, and angles-of-attack, a , at fixed Reynolds number. Re=2.xl&\ and thickness ratio, 
V;=59E are used to correlate the aerodynamic performance, measured by the power index, Q J/ VCy, which 
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appears explicitly in steady flight required-power equation. Aspect ratio and maximum camber form the 
input vector, p and C]‘ z / C D forms the target output vector, a, as shown below. 

AR 

\ y c 

where R is the number of input vectors of the training data. 

For the 3-D wing case, the maximum camber varies between 0.0 and 0.1 and the aspect ratio 
varies between one and Five. Three different training data sets are used out of the available data as shown in 
Table A 1 . Table A2 summarizes the test data sets used for prediction for this case. A simulation, referred in 
these tables, consists of two input variables: AR and y c and the output variable: Cl 3/2 /Cd- 

4. ASSESSMENT OF DATA PROCESSING AND OPTIMIZATION CAPABILITIES 

Of all the cases considered in this article, the impingement injector element, swirl co-axial injector 
element, 2-stage supersonic turbine and turbulent flow diffuser help understand the effectiveness of using 
polynomial based RSM. Shear co-axial injector element, 1 -stage turbine and two-dimensional wing model 
are used to carry out a comparative study between RSM and NN. The size of the data set used varies trom 
modest to large (from 9 to 2235 data). 

In the following, we synthesize the studies of Papilaet al. [49] & [50], Madsen et al.[33], Shyy et a!. 
[57], Tucker et al. [64] & [65] and Vaidyanathan et al.[71J. We First review the data processing capabilities 
then evaluate the performance of the optimization techniques. For both NN and polynomials, one needs to 
First decide most appropriate constructions for a given data set. For the NN, the choices are usually (l) the 
number of neurons, and (2) the error goals. Furthermore, the spread constant (for RBNN) and the number 
of hidden layers (for BPNN) can be specified. In this article, the BPNN and RBNN will be limited to the 2- 
layer form. 


= Cl 2 /C d L 


(30) 
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4.1 SHEAR CO-AXIAL INJECTOR 


4.1.1 Polynomial Fits 

According to the injector model developed by Calhoon et al. [3], injector performance, as 
measured by ERE depends only on the velocity ratio, V/V 0 , and combustion chamber length, L ctjm b . 
Examination of the original data set in Table A3-A6 indicates 15 distinct design points for ERE. Since 
chamber wall heat flux is dependent only on velocity ratio, V/V ( „ and oxidizer to fuel ratio, O/F, there are 9 
distinct design points for Q. The design space for this effort is depicted in Figure 14. For ERE , the 5 distinct 
chamber lengths offer the potential for a fourth-order polynomial fit in Lcombi while the three different 
velocity ratios limit the fit in V/V 0 to second order. Quadratic and cubic response surfaces for both ERE and 
Q have been generated for evaluation. These response surfaces represent reduced models accomplished by 
term elimination from the full surface using t-statistics as described earlier. The above-noted limitations on 
the data cause the cubic surfaces to be third order in L cumb only. 

Based on the adjusted RMS error, Vaidyanathan et al.[71] have concluded that the cubic fit is 
more accurate than the quadratic fit for ERE. The adjusted RMS error for the quadratic and cubic response 
surfaces of ERE are 0.211 and 0.083, respectively. The cubic fit, by this measure, is superior for ERE. 
However, the error is almost identical in the case of Q for both the quadratic (0.039) and cubic (0.040) 
surfaces, perhaps due to the very small number of points available for the curve fit. The additional terms in 
the cubic fit relative to the quadratic fit do not improve the mapping of the response surface for Q. 


4.1.2 Construction of RBNN 

In the case of the injector design there are two objectives, namely ERE and Q. Figure 15 gives the 
variation of crfor the network designed with solverbe for these objective functions. In case of solverb the 
error goal during training also defines the accuracy of the network. An objective of fitting a numerical 
model is to remove the noise associated with the data. A model, which maps exactly as solverbe does, will 
not eliminate the noise, w'hereas solverb will. Figure 15 shows that for low values of spread constant, the 
NN has a poor performance. As the spread constant increases a asymptotically decreases. However, as 
demonstrated by Figure 1 5a the performance of the network can deteriorate for higher values of the spread 
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constant . The region with a large variation in <ri s highly unreliable because this indicates a high sensitivity 
of the model to a small variation of spread constant and possibly the test data, in this region. Hence the 
desirable spread constant is selected from the region where the performance of the network is relatively 
consistent. 

Figure 16 gives the variation of <7 for the network design with solverb for the objective functions 
of ERE and Q. It also shows the influence of error goal on the network. Generally if a network maps the 
training data accurately it can be expected to perform efficiently with the test data. However, accurately 
mapping noisy data may result in poor prediction capabilities for the network. The variation in the 
performance is not significant except for the ERE and Q network (Figure 16), where the poor performance 
of the network at high values of spread constant improves for a larger error goal. This may indicate the 
presence of noise in the data for ERE , which solverb is able to eliminate with an appropriate error goal. 
Figure 17 shows variations in number of epochs and erwith the variation of error goal for a given spread 
constant when RBNN is designed with solverb. The number of neurons in the network is one more than the 
number of epochs. One expects that as the error goal increases the number ot epochs becomes smaller and 
the network performs less accurately as in Figure 1 7a and b. 

When choosing an appropriate network the above-mentioned features have to be considered. The 
performance of the constructed NN is best judged by comparing the prediction error as given in Eqn. (8) for 
different networks. Using solverbe , networks are designed with varying spread constants and the one that 
yields the smallest error is selected. When solverb is used, networks are designed for different spread 
constants and error goals. The network that gives the smallest error for the test data is used. The details of 
the networks selected are discussed in later sections. 


4.1.3 Evaluation of Polynomial and MV for Data Processing 

The pol>nomiai and NN-based RSM arc constructed using the training data. The test data is then 
employed to select the best polynomial or NN. Specifically in polynomial-based RSM, the difference 
between the polynomial and the training data, as given by cr a , is normally used to judge the performance of 
the fit. The additional use of the test data helps to evaluate the performance of different polynomials over 


34 



design points not used during the training phase. This gives a complementary insight into the quality ot the 
polynomial model over the design space. For both the rocket engine components, different polynomials 
were tried. Table 4 and 5 compare the performance of different polynomials used to represent the two 
objective functions of the injector case, ERE and Q. Starting with all the possible cubic terms in the model, 
revised models are generated by removing and adding terms. Similar kind of analysis is also done for the 
turbine case. The best polynomial is selected based on a combined evaluation between cr u and a 

For the NN, the test data helps evaluate the accuracy of networks with varying neurons in BPNN 
and varying spread constant in RBNN Thus the test data are part of the evaluation process to help select 
the final NN. Based on the RSM or NN model, a search for optimum design is carried out using a standard, 
gradient-based optimization algorithm over the response surfaces represented by the polynomials and 
trained neural networks. 

A reduced quadratic and an incomplete cubic response surfaces are used for the two objective 
functions. The first model in Table 4 and the sixth model in Table 5 are the selected cubic models tor ERE 
and Q< respectively. There is no noticeable improvement amongst the remaining cubic models for ERE. For 
2, the selected model is the best in terms of a x although there are other models with identical values of a 
The radial basis networks designed w r ith solverbe are the largest with 15 neurons in the hidden 
layer for ERE network and nine neurons for the Q network. Solverb designs a network for ERE with 14 
neurons in the hidden layer and a network for Q with eight neurons. Compared to RBNN, BPNN has fewer 
neurons, the number of neurons in the hidden layer are eight and four for the ERE and Q networks, 
respectively. Details of the networks used are listed in Table 6. The spread constant used for RBNN and 
the error goal of the training data is also given in this table. The spread constant values are selected from 
the region where the performance of the network is consistent with the variation of spread constant (Figure 
15 and Figure 16). The error goal, in the case of solverb , is selected based on the network with the best 

performance for the ideal spread constant. (Figure 17). 

The error in predicting the values of the objective function by different schemes is given in Table 

7. Several observations can be readily made. 

(1) Both NNs performs better than the RSM for this data set. 

(2) Both solverbe and solverb are of comparable performance. 
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(3) The BPNN helps generate smaller networks and hence performs at par in comparison to RBNN. 

(4) The cubic polynomial is more accurate than the quadratic one. 

The various models generated are compared with test data shown in Figure 18 and Figure 19. The 
curves representing the NN predictions are closer to the data obtained from the injector model than the RS 
thereby demonstrating that NN models are able to predict better than the RS. BPNN pertorms as well as 
RBNN but tends to be flat. Due to its lower order, the quadratic polynomial is tlat. The cubic polynomial is 
able to perform better than quadratic. 


4.1.4 Polynomial-based RSM for Design Optimization 

This case study is used to perform a complete comparative study between polynomial and NN- 
based RSM. The comparison is carried out in three ways. Firstly, the predictive capabilities of the different 
models are compared. Secondly, NN is used to increase the population of the design space, which is then 
used for mapping by polynomial -based RSM. Thirdly, polynomials and NN are used individually to 
represent the design space and help in the optimization of the design. 

An optimization was done for three different ranges of the independent variables using the 
quadratic fit. The three cases analyzed differ only in the constraints implemented on the design parameters. 
The constraints are 

Case I: 4 <0/F <6, 4 <V/V a <6. L comb <7. 

Case 2: 4 <0/F <6, 5 <V/V () <7, L (umb <7, 

Case 3: 4 <0/F <6, 6 <V/V a <8, L i0fnb <7. 

The optimization is repeated using the cubic fits. The combinations of weighting factors for ERE , 
5 , and for Q, t , are selected as (1,10), (1,1) and (10,1) for these three cases. The optimum has been 
evaluated and tabulated for each case, as detailed in Table 8 -Table 10. In this effort, injector element 
optimization means maximizing the performance, while minimizing heat tlux and chamber length. The 
optimum value for V/V t> obtained on the cubic response surface is quite different than that found on the 
quadratic surface for some cases (these particular cases are noted in bold in Table 8-10). Also, for selected 
cases where there are discrepancies between the quadratic and cubic results, the exact values from the 
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injector model have been included in parentheses in the tables for comparison. In these cases, the cubic fit 
more closely matches the exact data than does the quadratic fit. Sample results for ERE plotted in Figure 
20a clearly show the data is better fit by the cubic surface for the case shown. Figure 20b shows that the 
response surface predicted by cubic fit for Q has a noticeable dip that is completely missed by the quadratic 
fit. This discrepancy results in the optimum tor the cubic fit being considerably lower than that for the 
quadratic surface. The prediction from cubic fit agrees well with the exact data, which also has a dip for 
this specific case. 

The injector model was also used to produce additional design points to assess the capability of 
the different response surfaces to match the exact data. In Figure 21a and Figure 22a, the actual data 
obtained from the injector model for all the design points has been shown. The cubic and quadratic 
response surfaces obtained based on the original data is also shown. The RMS error lor predicting the new 
ERE data is 0.270 and 0.205 for the quadratic and cubic surfaces, respectively. For Q , it is 0.025 and 0.016 
for the quadratic and cubic surfaces, respectively. Again, the performance of the cubic surface is superior to 
that of the quadratic surface. 


4.1.5 Radial Basis Neural Networks (RBNN) 

Radial Basis Neural Networks are trained by both Solverbe and Solve rb for each injector design 
response, ERE and Q using the original data set of 45 design points. Solverbe trained the network for ERE 
with an error to the order of 1 0' 13 . The network trained by Solverbe for Q has an error on the order ot 10 l6 . 
Both networks represent the respective design spaces essentially exactly. Solverb , with an error goal of 
0.001, trained networks for both responses to represent the original data set adequately. Since the size of 
the data set considered for training the network is fairly modest, the number of neurons generated by 
solverbe is also small. Solverb would have been suited better for a larger data set where a reduction in the 
number of neurons might have appreciably reduced the computation time. The networks trained using 
Solverbe have been used for this article. The ability of the RBNN to fit the design data and to generate 
additional data for constructing a more accurate response surface is discussed in the follow-ing sections. 


37 



(i) Comparison between Solverbe and Solverb 

Since Solverbe trains with the same number of neurons (45 in this case) as data points, as seen 
above, it fits the training data set with negligible error. However, it can also create erratic behavior since it 
makes no attempt to filter noise generated by excess neurons in the network. Solverb , on the other hand, 
tends to reduce the potential for noise by controlling the number of neurons in the network. Table 1 1 shows 
that in the present article, for the spread constant value of 1.00, Solverb performs slightly better than 
Solverbe based on the nominal error measure. However, when judged by the level of errors associated, both 
RBNNs are satisfactory from a practical standpoint. As expected, Solverb uses fewer neurons than 
Solverbe ; in this case three less. It should be noted that, as investigated in detail by Papila et al. [49], the 
relative performance between Solverb and Solverbe is case dependent. 

(ii) Comparison of RBNN Predictions with Polynomial Based RSM 

Figure 21b and Figure 22b show that the RBNN trained by Solverbe is able to more accurately 
generate additional design data than either quadratic or cubic polynomial (shown for comparison in Figure 
21a and Figure 22a). In Figure 21a, the ERE surface trained with the original data set is shown. The 10 
extra design points calculated with the injector model for V/V 0 of 5.00 and 7.00 are shown. The ability of 
the RBNN to accurately generate new design data can be seen by comparing the fit for ERE in Figure 21b 
to that for the polynomials in Figure 2 1 a. RBNN trains the network with more flexibility and learns the data 
trend, whereas polynomials provide only an approximate fit on the given data. Regarding the RMS error, a, 
for ERE , it is 0 1 52 for RBNN predictions as compared to the values of 0.270 and 0.205 for quadratic and 
cubic surfaces, respectively. The four extra design points generated for Q t also at V/V a of 5.00 and 7.00 are 
shown compared to the polynomial surface in Figure 22b and compared to the RBNN surface in Figure 
22b. The RMS error in the case of Q is 0.022 for RBNN as compared to 0.025 and 0.016 for quadratic and 
cubic surfaces respectively. Here the performance of the RBNN is better than the quadratic but slightly 
poorer than the cubic fit. Examination of Table 1 1 indicates it may be possible that using Solverb with a 
spread constant of I 05 could further reduce the RMS for Q. However, the errors for Q are low enough that 
further reduction may not be practical. 
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4.1.6 RBNN-Enhanced Polynomial Based Response Surface 

Additional design points generated by the RBNN are added to the original data set to form the 
enhanced data set. This enhanced data set is used for further analysis to evaluate the performance of the 
RSM with the larger number of design points. The enhanced data set for ERE has 15 points from the 
injector model and 10 from the RBNN, for a total of 25 points. The enhanced data set for Q has 9 points 
from the injector model and 4 from the RBNN, for a total of 13 points. The entire optimization analysis was 
redone with the enhanced data set. On this enhanced data set, the full quadratic response surface seems 
already appropriately constructed and invoking the statistical analysis generates no reduced model. With 
the added data in the enhanced data set, it is now possible to obtain a fit for ERE that is 4 th order in V/V () 
and 4 th order in L comb . Q can now be fit with a cubic in V/V a and a quadratic in O/F. This is now possible 
since a combination of 3 different values of O/F, 5 different values of V/V a and 5 different values ot L com i, 
are available. 

(i) Comparison of Fits with the Original Response Surfaces 

Comparison of the enhanced response surfaces with the original response surfaces indicates that 
the extra data produced with the RBNN generally improves the quality of the curve fit. The adjusted RMS 
error for ERE on the original set is 0.211 and 0.083 for quadratic and cubic fits, respectively. On the 
enhanced data set, it is 0. 179 and 0 100 for the quadratic and cubic fits, respectively. The slight increase in 
the error in the case of the cubic fit may be due to noise related to the over-sensitivity ot the polynomial. 
However, this phenomenon may reflect the fact that the level of the RMS is low- enough in either case so 
that no further improvement is accomplished. The adjusted RMS error for Q with the original set is 0.039 
and 0.040 for the quadratic and cubic fits, respectively. On the enhanced set it was 0.027 and 0.026 for the 
quadratic and cubic, respectively. With the exception of the cubic fit for ERE , the fits from the enhanced 
surface are improved over those from the original surface. Also, when optimum design points are 
examined, there is less difference between the quadratic and cubic fits on the enhanced surfaces than there 
is on the original surfaces. 
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(ii) Comparison of Optimal Design Points 

The analysis for the three cases of optimization over the same three ranges of independent 
variables has been re-done. The results of the optimization on surfaces generated from the enhanced data 
set are tabulated in Table 12-Table 14. The predicted optimal design points using cubic and quadratic tits 
are generally close to each other. They are closer to each other on the enhanced data set than on the 
surfaces generated using the original data set. One case where the cubic and quadratic optimum points are 
somewhat different is analyzed further. The results shown in Figure 23 confirm the optimum value of 
velocity ratio on the quadratic fit to be low-er than the cubic fit in this case. Given the weightings of 1.0 for 
ERE and 10.0 for Q , the optimizer has selects the minimum of Q for both fits. Since the curves exhibit 
different minimum points, the weightings force the selection of different optimum points. As already 
discussed, for the polynomial fits on the RBNN-enhanced data sets, the errors of both quadratic and cubic 
polynomials are more comparable than in the original analysis. At the upper limit ot the design space tor 
combustor length, the ERE curves tend to flatten out. This causes some difficulty in locating the optimum 
and may cause more noticeable differences between the different polynomials. However, different optimal 
designs selected by different polynomials under such a circumstance are not important since these yield 
very similar injector performance 

The optimum solution obtained from various schemes is shown in Table 15 and Figure 24 and 
Figure 25. The aim is to maximize ERE and minimize Q. The trend of the objective functions in the design 
space is monotonic and hence every model is able to select identical optimum design for the given 
constraints. The flatness of the quadratic polynomial results in less accurate values of the objective function 
for the optimum design. The cubic polynomial, while more flexible than quadratic, is not consistently better 
in predicting the optimal design point. For example, a V/V 0 constraint of 4, the quadratic polynomial is 
more accurate but for higher values of V/V„ the cubic polynomial is more accurate. In contrast, the NN 
models are able to perform well. Since the optimum design happens to be the same as one of the training 
points, solver be is able to predict the values of the objective function accurately. Solverb performs equally 
well, thereby showing the capability of performance with fewer neurons. Performance of BPNN is not as 
satisfactory as suggested in Table 7. For lower constraints of V/V ( „ it performs poorly but for higher values 
of V/V n it is good. This may be due to the selection of fewer neurons in the hidden layers of the networks. 
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Overall, it is still better than the polynomial -based RSM and demonstrates the flexibility of NN over 
polynomials. 

As stated by Papila et al. [49], when it comes to choosing between NN and polynomials, 
polynomials are easy to compute. The number of coefficients might be numerous but the linearity of the 
system expedites the process of coefficient evaluations. This is also the reason RBNN train fast. On the 
other hand, the weights of BPNN are evaluated through a nonlinear optimization, which slows the training 
process. Of all the NN presented here, the one designed w r ith the help of solverbe is the fastest to train since 
the values of the weights are set to values of the input dependent variables. Solverb trains with the addition 
of one neuron at a time with weights similar to the input and hence is slower. 

4.2 IMPINGEMENT INJECTOR ELEMENT 
4.2.1 Polynomial Fits 

In Tucker et al. [64], method i uses the polynomial based RSM to find optimal values of ERE , Q w , 
Q inj , W re t and C re i for acceptable values of AP (J , A P Jt and a The approach of RSM is to perform a 

series of experiments, or numerical analyses, for a prescribed set of design points, and to construct a 
response surface of the measured quantity over the design space. In the present context, the live responses 
of interest are ERE , Q inp W rei and C'V The design space consists of the set of relevant design variables 
AP ()> A P h L, omh and a 

(i) Individual Polynomial Models 

When the JMP software is used to analyze the 163 design points, five individual full response surfaces for 
the variables in the design space are approximated by quadratic polynomials that contain 15 terms each. 
Using the t-statistics approach noted above and detailed in Tucker et al. [64], unnecessary terms in each 
equation can be eliminated to give the reduced quadratic surfaces A survey of the reduced response 
surfaces indicates that the equations reflect the functionality used to construct the models for the dependent 
variables. 
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(ii) Joint Response Surfaces 


In the current article, it is desirable to attempt to maximize ERE and while simultaneously 
minimizing Q w , Q mj , W rei and C re {. Therefore composite response surface for the present case is given by: 


/ \t/5 

^ ( d ERE d Q ^ d w t d c ^ j 


(31) 


4.2.2 Optimization Results & Discussion 

Three sets of results are presented below to demonstrate the capability of method i for the current 
injector design. These three examples illustrate the effect of each variable on the optimum design, the 
trade-offs between life and performance issues, and the effect on the design of extracting the last increment 
of performance. 

(i) Effect of Each Variable on the Design Using Original Constraints & Equal Weights 

The results in this section were obtained by building the joint response surface with the addition of 
one dependent variable at a time. The results are shown in Table 16. Since current non-optimizer based 
design methods yield high-performing injector elements, simply maximizing the ERE is not a challenge. 
Accordingly, the initial results (Case 1 ) are obtained with a joint ERE and Q w response surface. The results 
in Case 2 have the impingement height added. Case 3 adds the relative chamber weight and the relative cost 
is added in Case 4. All results are obtained using the original independent variable constraints and all 
dependent variables have equal weights of one. The results for Case l show that ERE is at its maximum 
and Q w is very near its minimum desirability limit. Minimizing Q w requires a small AP } relative to AP a as 
evidenced by the values of lOOpsi and 183psi, respectively. Maximum ERE values are found at the longest 
chamber length, U rr „„/>=8inches. Even with the relatively high value of 183psi for AP 0 and low value of AP f 
of lOOpsi, ERE is maximized to 99.9% with an impingement half-angle of 33. 1°. 

Addition of the impingement height to Case 2 to model the injector face heat flux, Q m} , forces a 
lower to increase H impirtf{e and decrease Q inJ . This decrease in the radial component of the fuel momentum 
has an adverse affect on ERE. This effect is mitigated to a degree by increasing the AP S by 32psi to 132psi. 



ERE is still reduced by 1.6%. Also, the increase in AP causes increased penetration of the fuel jet, which 
results in a slightly higher Q w . 

Case 3 adds the relative combustor weight to the list of dependent variables modeled. Since W re( is 
only a function of L comby minimizing W rfi shortens the combustor length from 8 to 6.6 inches. The shorter 
L CO mb tends to low'er ERE. This effect is offset to a large degree by increases in ziP, and ca both of which 
increase the radial component of the fuel momentum. The increase in AP, also causes a slight increase in 
Q w . The increase in causes a significant decrease in H impinfiey which increases the injector face heat flux. 
Finally, the relative cost of the injector is added in Case 5. Since C re ( is only a function of propellant 
pressure drops, both AP 0 and AP f are driven to their respective minimum values. This and a slight increase 
in a allow ERE to be maintained at 98%, even with a slight decrease in L como . The largest effect of this 
fairly dramatic decrease in propellant pressure drops is on Q w . Even though the values for AP 0 and AP f fell, 
AP, increased relative to AP„ causing Q w to increase by almost 9%. Impingement height and relative 
combustor weight are essentially unchanged. 

Although several of the variables included in this exercise are qualitative, an important conclusion 
can still be drawn. The sequential addition of dependent variables to an existing design results in changes to 
both the independent and dependent variables in the existing design. The direction and magnitude of these 
changes depends on the sensitivity of the variables, but the changes may well be significant. The design in 
Case 4 is quite different that the one in Case l. Consideration of a larger design space results in a different 
design — the sooner the additional variables are considered, the more robust the final design will be. 

(ii) Emphasis on Life & Performance Issues Using Original Constraints & Unequal Weights 

The purpose of this section is to illustrate the effect of emphasizing certain design criterion on the 
optimization process. Method i allows this emphasis via the weights applied to the desirability functions in 
the joint response surface. The results shown in Table 17 facilitate the illustration. The Case 1 (baseline) 
results are repeated from Case 4 in this table where the entire design space is considered with the original 
constraints and equal weights for the dependent variables. The results in the Case 2 column are obtained by 
emphasizing the minimization of the wall and injector face heat fluxes. Desirability functions for both of 
these variables are given a weight of five. Since lower heat fluxes tend to increase component life, 
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weighting these two variables is equivalent to emphasizing a life-type issue in the design. As expected, a is 
decreased to increase H (mpinge , thus decreasing Q :n; . Since the fuel pressure drop is already at the minimum, 
the oxidizer pressure drop is increased by 58% to decrease Q w . Both of these changes tend to decrease ERE 
While ERE does decrease, the effect is somewhat mitigated by an increase in L comb . The increases in Uomb 
and AP 0 cause increases in W rel and C reh respectively. The emphasis on life extracts the expected penalty on 
performance. Additionally, for the current model, there are also weight and cost penalties. 

The results for Case 3 are obtained by emphasizing maximization of£7?£and minimization of W re! 
with desirability weightings of five. Increased weighting for these two variables is equivalent to 
emphasizing a thrust to weight goal for the injector/chamber. The relative chamber length is shortened to 
lower W rel . ERE is maximized by increasing the radial momentum of the fuel jet. Both APj and a are 
increased to accomplish ERE maximization. As noted earlier, increasing dlP, and a lead to increased wall 
and injector heat fluxes, respectively. Table 17 indicates that to be the case here. For this case, emphasis on 
thrust and weight tend to have an adverse affect on both Q w and Q inj . Relative cost, tor the current model, is 
not significantly affected. 

(iii) Extraction of Last Performance & Weight Increments (Modified Constraints & Unequal 
Weights) 

Here, the high marginal cost of realizing the last increment of thrust to weight is shown. This 
section illustrates the capability to modify the constraints on the independent variables and use unequal 
weights on the dependent variables at the same time. The results for Case 3 in Table 17 are carried over to 
Case 1 in Table 18 as the baseline for this example. Here the original constraints are used but increased 
weights have been applied to emphasize ERE and VV rW . Cases 2 and 3 modify the constraints on the 
propellant pressure drops, raising the minimum pressure drop from lOOpsi to ISOpsi. For Case 2, both AP i} 
and zlP/ are now at the minimum level for the modified constraints. U- 0 mi> is increased slightly to maintain 
ERE. The decrease of AP, relative to AP rt causes a decrease in Q w . The slightly higher-pressure drops also 
cause Crei to increase somewhat. Other variables are not changed appreciably. 

For Case 3, ERE and VV rW are further emphasized by increasing their desirability weights to 10 
while decreasing the other weights to 0. 1. L itmh is shortened to respond to the increased emphasis on weight 
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minimization. Maintaining the high level of ERE requires large increases in dP f and ato increase the radial 
component of the fuel jet momentum. The increase in APf causes over-penetration of the fuel jet, which 
results in an increase in wall heat flux. The large increase in a yields the expected decrease in 
which increases the injector face heat flux. The additional emphasis on ERE and C ret yields essentially no 
increase in ERE in this range of AP, although a small weight savings is seen. These marginal improvements 
are offset by fairly large increases in C r <i and Q inj . 

4.3 SWIRL CO-AXIAL INJECTOR ELEMENT 

Two sets of results are presented below to demonstrate the capability and flexibility of method i 
for the current injector design. These examples illustrate the effect of each variable on the optimum design 
and the trade-offs between life and performance issues. 


4.3.1 Effect Of Each Variable On Element Design 

The results in this section were obtained by building the joint response surface with the addition of 
one dependent variable at a time. The results are shown in Table 19. Case 1 seeks the maximum 
performance without regard to the effect on the other dependent variables. ERE is a fairly strong function 
of L romt , — longer chamber lengths allow more residence time for the propellant to mix and burn. The effect 
of <9 on ERE is strongest at low values of 0 ERE increases with increasing 0 until about <9=80° and then 
fall off slightly due to the competing influences noted earlier. These competing influences also cause the 
effect of both pressure drops on ERE to be somewhat flat, although since AP 0 affects more variables, its 
influence is slightly stronger. Maximum performance is found at high values of AP a , 0 , and Lc omb and at 
low values of J P ; . This trend is consistent with other works for similar injector elements. The predicted 
optimal value of 98.5 is indeed the highest predicted by this model. 

The objective of Case 2 is to simultaneously maximize ERE and the minimize Q w . Table 19 shows 
that the exact same design point was chosen as for Case 1. Usually, the design, which yields the maximum 
ERE, also produces a high wall heat flux. That is not the case here; this issue has already been noted. The 
minimum Q„ is found in the region of high AP t) and low J P f . In this area, Q w is almost independent ot 0. 
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Hence, the minimum Q w can still be found for a high value of <9 required to maximize ERE. It should be 
noted that in the low A P<„ high AP f region, Q+ is a function of & Here, as O is increased, Q w increases 
since the larger swirl angle forces d Q to increase and thus decrease V 0 . In the Calhoon et al. [3] model, this 
reduction in G0 2 momentum causes an increase in Q„. 

The requirement to minimize Q lftJ is added in Case 3. In order to minimize Q tnJ , the swirl angle is 
decreased from 81° to 37°, thus reducing the injector face heat flux by approximately a factor of 3. This 
decrease in <9 also lowers ERE which forces use of a longer chamber to offset some of the loss. Still, ERE 
is reduced by over one percent. 

Case 4 considers the desire to minimize the chamber weight, W reh in addition to maximizing ERE 
and minimizing Q w and Q inJ . Since W re , depends only on L C()mh , the chamber length is shortened by over half. 
The weight is reduced, but so is ERE . To mitigate the adverse effect on ERE , <9 is increased by almost 10°, 
simultaneously increasing Q :n; . ERE drops again by over a percent, while Q w remains constant. 

Finally, minimizing the injector cost. C rei , is added in Case 5. Decreasing each pressure drop 
approximately a factor of 2 lowers the relative injector cost. Decreasing ^ results in a larger fuel gap and 
decreasing AP a allows for a larger swirl slot. These factors combine to lower the cost by almost 10 %. 

Although several of the variables included in this exercise are qualitative, an important conclusion 
can still be drawn. The sequential addition of dependent variables to an existing design results in changes to 
independent and dependent variables in the existing design. The direction and magnitude ot these changes 
depends on the sensitivity of the variables, but the changes may well be significant. The design in Case 5 is 
quite different that the one in Case 1. Consideration of a larger design space results in a different design — 
the sooner the additional variables are considered, the more robust the final design. 


4.3.2 Emphasis on Life and Performance Issues 

Method i allows this emphasis via the weights applied to the desirability functions in the joint 
response surface. The set of results shown in Table 20 facilitate the illustration. The baseline results Table 
20 (repeated from Case 5 in Table 19) consider the entire design space using the original constraints and 
equal weights for the dependent variables. The results are obtained by emphasizing the minimization of the 
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wall and injector face heat fluxes for Case l. Desirability functions for both of these variables are given 
increased weights (5 and 10, respectively). Since lower heat fluxes tend to increase component life, 
weighting these two variables is equivalent to emphasizing a life-type issue in the design. Since Q w is 
already at its minimum value, it remains fixed. As expected, 6>is decreased, which decreases the value ot 
Q inj by almost 35%. The lower value of <9 also produces a lower ERE. Both propellant pressure drops and 
the combustor length are increased to mitigate the drop in ERE. The increases in L comb and AP, cause 
increases in W re{ and C reh respectively. The emphasis on life extracts the expected penalty on performance. 
Additionally, for the current model, there are also slight weight and cost penalties. 

The results for Case 2 are obtained by emphasizing maximization of ERE and minimization of W re{ 
with desirability weightings of 10 and 5, respectively. Increased weighting for these two variables is 
equivalent to emphasizing a thrust to weight goal for the injector/chamber. The relative chamber length is 
shortened to slightly lower W rel . ERE is maximized by increasing the G0 2 swirl angle by a factor ot almost 
2.5 and also increasing AP t by over 35 %. The value of ERE rises by over one percent. As noted earlier, 
increasing <9 leads to increased injector heat flux. For this case, emphasis on thrust and weight tends to 
have an adverse affect on Q inj . Relative cost, for the current model, is also increased significantly. 
Performance and weight trends for the swirl and impinging injector elements are shown in Figure 26. 
Figure 27 shows the heat flux and cost trends for the swirl and impinging injector elements. 

4.4 SUPERSONIC TURBINE FOR REUSABLE LAUNCH VEHICLES 

4.4.1 Polynomial-Based RSM Results for 1-, 2- and 3-Stage Turbines 

There are 28-unknown coefficients needed for constructing the 2 nd -order response surface for the 
single-stage case, 78 for the 2-stage and 136 for the 3-stage case. Different starting points are tried to avoid 
local maximum and the optimum values of t}, W and Apay with the corresponding design parameters are 
determined. The results shown are comparable with the corresponding Meanline runs with the highest error 
of 5% for Apay for single-stage turbine. The percentage error is increased to 13.5% for Apay for 2-stage 
turbine and to 14.6 % for the 3-stage turbine for Apay indicating that the accuracies of the response surfaces 
constructed are poor for the 2- and 3-stage. 
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Papila et al. [50] have reduced the size of the parameter space by 80% in each coordinates, based 
on the optimal design selected in the original design space, to improve the accuracy of the response 
surfaces for these cases. The intention is to improve the fidelity of the response surface. With these refined 
designed spaces, substantial improvement of the response surface fit accuracy is observed for both cases by 
Papila et al. [50]. 

Based on the results obtained, the following observations can be made: 

(1) To ascertain required predictive capability of the RSM, a two-level domain refinement strategy 
has been adopted by Papila et al. [50]. The accuracy of the predicted optimal design points based on this 
approach is shown to be satisfactory. 

(2) According to the results obtained for Apay-based optimization, the 2-stage turbine gives the best 
Apay result. 

(3) As the number of stages increases, it is observed that efficiency improves w-hile the weight 
increases also but the improvement in efficiency can’t compensate the penalty from higher weight. 

(4) As shown in Figure 28, the mean diameter, speed, and the exit blade area exhibit distinct trends. 
Specifically, the diameter decreases, speed increases, and annulus area decreases w r ith increasing number of 
stages. It is interesting to observe that none of these design parameters are toward the limiting values, 
indicating that the optimal designs result from compromises between competing parametric trends. For 
such cases, a formal optimizer such as the present response surface method is very useful . 


Table 21 gives a summary of the optimization results for 1-, 2- and 3-stage turbines for Apay- 
based optimization. 


4.4.2 Higher Order Polynomials and NN-Based RSM for Single-Stage Turbine 

The generation of polynomial-based RS model and the training of the NN are done with 76 design 
points of the single-stage turbine. The analysis was initially done without the constraints and then with the 
constraints on (AN) 2 and V pucft . 

A quadratic polynomial model was initially generated. Then, cubic terms were included. Cubic terms that 
are products of three different variables were included because of the number of data available and the 
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number of levels being three. The trend of the design data also suggests the presence of some of these 
terms. Therefore, the initial cubic equation has 45 terms. Reduced third order polynomial model for r) and 
W were selected based on the relative performances of different polynomials obtained by removing terms 
from the initial cubic equation based on t-statistics. The cubic equation was selected based on the evaluated 
value of cr a and a Table 22 suggests that the reduced cubic polynomial is better than the quadratic 
polynomial since <r a is better for the former. The value of c is comparable. 

When constructing the NN-based response surface, the design parameters of the NN should be 
selected carefully since the selection of the design parameters determines the learning characteristics of the 
NN. For the single-stage supersonic turbine case, the variation of a with respect to the only design 
parameter of solverbe network, spread constant, is plotted in Figure 29 for both objective functions of rj 
and W. Figure 30 shows that for low values of spread constant, the NN has a poor performance. As the 
spread constant increases <7 asymptotically decreases. Therefore, the appropriate spread constant is 
selected from the region where the performance of the network is relatively consistent. Figure 3 1 shows the 
influence of error goal on the network performance. Unlike the case of injector (Figure 17), a more 
stringent error goal for the training data does not necessarily result in better predictive capability against 
the test data for the single-stage turbine. 

The networks designed with solve rb have 37 and 75 neurons tor rj and W, respectively in the 
hidden layer, while those designed with solverbe has 76 neurons each. The BPNN uses significantly less 
number of neurons by generating networks with five and 60 neurons for r) and W y respectively, in a single 
hidden layer. The NN architectures chosen are listed in Table 23. 

The accuracy of the various models is tested with the 18 additional available data and the error is 
shown in Table 24. Solverbe has a poor prediction for rj, which might be due to over fitting, but performs 
well for W. Solverb is most consistent among all methods evaluated. 

The optimum solutions subjected to the constraints, of (AN) 2 limited to less than 1.132 
(normalized with baseline value) and V puch is limited to less than 1.148 (normalized with baseline value), 
are presented in Table 25. Since (AN) 2 is proportional to the product of square of RPM and A ann , and V pt(ch is 
proportional to D times RPM , no NN or polynomial -based RSM is generated for them. By comparing the 
predicted optimal design by the various methods, one observes that solverbe and BPNN yield noticeably 
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larger error in rj and W y respectively. Solverb and the response surface are more consistent with both rj and 
W. Judged by the error in predicting Apay , it seems that the polynomial based RSM is most accurate. 
However, since the real goal is to maximize Apay , it is important to note that the actual value ot Apay lor 
the optimal design chosen by the RSM is the worst. 

From a design perspective, it is interesting to understand the impact of the constraints from A ann 
and V pitch on the optimal turbine parameters. Such an assessment is offered in Figure 32 and 33. As D. RPM 
and A ann decrease, rf, W , V pilc h, AN 2 and Apay decrease. C b and C v are almost constant over the design space 
and they do not have any noticeable effect on the objective functions and constraints. In the case of Cv, the 
BPNN shows a small perturbation for the analysis with the constraint. This might be due to the mapping ot 
some noise by BPNN. Otherwise it is unaffected by the inclusion of the constraints. The optimum stage 
reaction, K n is zero implying that the optimum design should be an impulse turbine. 


4.4,3 Orthogonal Arrays For 2-Stage Turbine 

Although the majority of the work is based on the face centered composite design approach 
(FCCD), orthogonal arrays (OA) are constructed by Papila et al. [50] to investigate the efficiency of 
orthogonal array designs in representing the design space tor 2-stage turbine. A set ol 249 design points is 
selected usins orthogonal arrays. Table 26 shows the comparison ol the quality ol the second-order 
response surfaces generated for r|, W and Apay by using 1990-data generated by face centered composite 
design and 249-data selected by orthogonal array method. 

The above table illustrates that the fidelity of the response surface generated for design space of 
249 data, based on orthogonal arrays, are comparable with that of 1990 data based on the face centered 
criterion. The response surface models are also assessed by using 78-test data to determine the predictive 
accuracy of these models. Table 27 presents that the testing adjusted rms-errors of response surfaces 
generated are 1.65% for >7 and 0.96% for W using 249-data, and 1.67% for /?and 1.21 9c for W using 1990- 
data. 

When these results are compared with the results of 1990-data and it is observed that the optimum 
7 , W and Apax are largely consistent However, it is also observed from Figure 34 which shows the 
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comparison of the design variables for optimization based on (Apay), some of the design variables are 
different even though optimum rj , W and Apay are consistent. This shows that there are multiple points in 
the design space, which yield comparable performance. Nevertheless, it remains true that the two-stage 
turbine is most suitable from a payload point of view. 


4.4.4 NN-Based RSM for 2-Stage Turbine 

In order to find the optimum RBNN design for the design of the two-stage turbine design, the 
effect of the spread constant (sc) on the network training error is determined. Figure 35 and Figure 36 show 
the variation of solverbe network error, a, with respect to spread constant for the NN designed for FCCD 
and OA data. The optimum spread constant is determined as 3.2 for 1990-training data (FCCD) from 
Figure 35 and 4.3 for 249-data (OA) from Figure 36. In spite of the fact that, the spread constants larger 
than 3 give reasonable training rms-errors (<7j, less than 0. 1% for all networks designed for refined space 
with 249-data as shown in Figure 37, 5c=4.3 value is used for these cases for consistency. 

After constructing the NN-based response surface, the NN model is tested by using 78-test data 
selected along the main diagonal of the design space to determine the predictive accuracy of these models. 
Table 28 presents that the prediction rms-errors (a) of response surfaces generated by for second order 
polynomial are 1.659c for t] and 0.9691- tor W using 249-data, and 1.67% tor rj and 1.21% tor W using 
1990-data. Table 28 also presents that the prediction rms-errors of response surfaces generated by solverbe 
RBNN are 1.36% for r} and 1.30% tor W , and 2.26 % for r} and 1.56% for W using 249-data. 

Figure 38 summarizes fitting/training and testing results of RBNN and polynomial -based Apay 
approximations for 2-stage turbine. The efficiency of the multi-level RSM approach can be observed by 
comparing the original and refined design space plots. From these plots, it is also possible to observe that 
more accurate training is possible with RBNN but testing or prediction accuracies of the RBNN and 
polynomial -based approximations are quite comparable. 


51 


4.5 TURBULENT PLANAR DIFFUSER 


4.5.1 Polynomial Fits 

Based on the D-optimal set of 35 design points selected, the 21 regressors of a full quadratic 
polynomial were fitted resulting in a moderate /? d 2 -v alue of 0:8 10. A backward elimination ot regressor 
terms subsequently led to the removal of five terms and an increase of R a 2 to 0.848. The lower values ot 
R a 2 , in comparison to the two-design-variable case, reflect the increased difficulties in obtaining a good fit 
when moving to higher-dimensional response surfaces. Data on the backward elimination steps are given in 
Table 29, which apart from R 2 and R a 2 holds the minimum t-statistic and the number of uncertain terms 
with |to| < 2.0 remaining in the model. From the t-statistics information, it appears that the backward 
elimination improved the accuracy of remaining terms. 

The next step performed was to investigate whether the 35 applied observations included outliers. 
A common (but not necessarily true) assumption, which enables the statistical treatment of observations, is 
that errors are independently and identically distributed according to a normal distribution with mean zero 
and variance. Thus, the distribution of response surface errors was plotted and compared to a normal 
distribution, with which it is expected to correspond well. From the histogram plot of the error distribution, 
see Figure 39, it did not seem that there are any outliers. Four arbitrary points away from sampling points 
were picked to test the prediction accuracy of the polynomial-based RSM. Table 30 compares CFD-results 
and polynomial approximations with and without backward elimination of terms. 

Again, the predictions of the response surface appear reliable, apart from at the last control point. 
This point is, however, well in the non-monotonic region, so that the approximation relies on an 
extrapolation, which was never intended. The reduced approximation model comes closer to the CFD- 
results for two out of the three meaningful test points. 


4.5.2 Numerical Noise 

While noisy data from laboratory experiments is a generally accepted fact, the presence of noise in 
numerical simulations seems much less recognized. Due to the complex numerical modeling techniques of 
CFD, the exact origins of noisy responses may be difficult to pinpoint, but factors such as turbulence 
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models, incomplete convergence, and the discretization itself are certainly influential. Here, the presence of 
numerical noise has been investigated. The problem of non- smooth or noisy objective functions has 
previously been addressed by Giunta et al. [17], who found RS approximations-based optimization to 
perform very robustly under such circumstances, especially when point selection is based on design of 
experiment techniques, such as D-optimal designs. 

Limitations of the software used were felt during the application of a wall shape parameterization 
in the investigation of noise. A B-spline curve with two free control points was used. Again, it was 
observed that the objective function oscillated due to numerical noise, but the amplitude was small. To 
make the noise more apparent, it was therefore necessary to refine the subdivision of the discretized line 
and reduce its length to 20% of the initial, so that the line spans from (0.3, 0.6) to (0.302,0.602). This 
yielded the noisy response patterns shown in Figure 40. The two topmost curves in this figure were 
determined using a relatively tight convergence criterion, and two different convection schemes - a 
standard first-order upwind differencing scheme (UDS) and a second-order upwind differencing scheme 
(SUDS). The use of different differencing schemes w-as carried out to estimate whether numerical diffusion 
does significantly dampen the generation of noise. As discussed in Madsen et al. [33], two different CFD 
codes w-ere adopted, and one seems less forgiving, in the sense that it predicts a stronger tendency for flow 
separation. This could possibly be explained by factors such as numerical diffusion, boundary treatments, 
and momentum interpolation methods adopted in the two codes. 

As expected, switching to a more dissipating differencing scheme (lower order accuracy) yields a 
smoother response. To further illustrate this issue, one more design line curve is shown in Figure 40, which 
arose from using a relatively loose, yet still reasonable, convergence criterion (using SUDS). The applied 
convergence criterion considers summed and normalized (by inlet flux) residuals over the entire mesh, with 
termination of computations once the maximum is below a certain small value e. The loose convergence 
criterion in Figure 40 was £= 10 \ whereas the tight tolerance was e = 10 5 . For comparison, a convergence 
limit of £ = t0 4 was applied in the CFD analyses used for response surface modeling. The overall 
conclusion is that the presence of some numerical noise in CFD-results is practically inevitable, although 
its magnitude depends on choice of code and modeling techniques. Here, a technique such as polynomial- 
based RSM can be effective in smoothing out the undesirable fluctuations. 
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4.5.3 Optimum Diffuser Designs 

In the optimum design using B-spline parameterization, both the monotonicity constraint and four 
out of five side constraints are active. As already mentioned, the response surface constructed to guarantee 
wall monotonicity becomes too restrictive. To compensate for this, a one-dimensional search in the 
direction of the steepest gradient was conducted starting at the optimum design point estimated by RSM. 

y = y* +aVF (32) 

This search is terminated as soon as designs turn non -monotonic, yielding a new optimum point at 
the edge of the true feasible domain and an increase in the optimum pressure recovery coefficient from 
0.7208 to 0.7235. Figure 41 compares the optimum wall contours determined by RSM using B-splines and 
polynomial shapes. The optimum B-spline shape compares well to the optimum polynomial one, so it is not 
surprising that there is no significant gain compared to this case. The largest differences in shape are found 
in the later part of the expansion, where the shape has less impact on the overall performance, as separation 
is small in either case. Thus, the close resemblance of optimum inlet shapes is reassuring in terms of the 
credibility of the optimization algorithm. A CFD-analysis of the five-design-variable optimum design 
yields a pressure recovery coefficient of 0.7193, a little below the predicted value, as in the two-design- 
vanable case. The improvement from the two design variable case (0.7185 to 0.7 193) indicates that there is 
not much potential for further gains. Furthermore, for comparison. Figure 41 also contains the 
corresponding wall contour determined using search optimization techniques. The optimum wall shape 
found by search optimization can be described as truly bell-shaped, without a "plateau” similar to the one 
found in the results of RSM-optimization. There appears to be a distinct difference in optimum shapes from 
the two different optimization approaches, which must be ascribed to the combination of optimization 
accentuating modeling differences and a relatively small scatter in diffuser performances. 

Figure 42 highlights the use of a response surface approximation for the optimum shape of a two- 
dimensional diffuser. As illustrated, within the fidelity of the analysis tool, there are often multiple design 
points that meet the design objectives. It is interesting to note that different diffuser shapes are found to 
yield essentially the same performance. The response surface model is ideally suitable for such situations. 
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4.6 LOW REYNOLDS NUMBER WING MODEL 


4.6.1 Polynomial Fits 

For the 3-D wing case, the response is the flight power index, C*' 2 IC 0 , and the design space 

consists of design variables maximum camber, y, , and wing aspect ratio, AR Quadratic, cubic and quartic 
order polynomials are tested for the best approximations for data sets containing 9, 15 and 25 simulated 
data points (See Table Al). The predicted RMS errors are calculated for each of the model and are shown 
in Table 31. As shown in this table, Model 4 gives the smallest predicted RMS for the cases involving 9 
and 15 simulated data points, whereas, Model 12 allows the smallest predicted RMS error for the case 
involving 25 simulated data points. 


4.6.2 Comparison of Radial-Basis and Back-propagation Networks 

The predictive accuracy of neural networks depends not only on the training data but also on the 
parameters used to define the network. The best values for these parameters cannot be determined by using 
only training data, because typically one can obtain very small errors for the training data with a wide range 
of these parameters. However, the performance of NN can be examined using test data. 

For the radial-basis network, one important issue is to investigate the magnitude of error in the test 
data to help to select the spread constant. For the back-propagation network, where cost of computation is 
an issue, the effect of number of neurons on the cost and accuracy should be checked. It was noticed that 
for the back-propagation network, using four neurons gave a good compromise of accuracy and cost. For 
the radial-basis network, it was found that the error and the number of iteration required for convergence 
are extremely sensitive to the value of spread constant. After extensive experimentation, the spread 
constant was chosen as 1. 175. 

For the 3-D wing case, both radial-basis NN and back-propagation networks are applied. In order 
to be able to make comparisons between these networks, the training time histories are summarized in 
Table 32 and Table 33. These tables show' that both are efficient in the training of 9-simulation, 15- 
simulation and 25-simulation training data sets. However, as the data size increases, the back-propagation 
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network exhibits a growth rate in terms of the number of epochs, indicating that it is more CPU time 
intensive for larger data sizes. As far as accuracy is concerned, both networks perform well exhibiting 
improved predictive capabilities as the number of training points increases from 9 to 25-simulation for y c 
interpolations ( Figure 43). For this case, both methods reproduced the original 9-simulation accurately but 
both failed to predict accurately the interpolation points at y c =0.0 125,0.025, 0.075 and 0.0875 with the 
RMS error of the test data of 1.68 for back-propagation network and 1.04 for radial-basis network ( Figure 
43a). Figure 43b shows that adding 6 new points at AR = 2 and 4 at v t =0, 0.05 and 0.1 (15-simulation 
training data set) does not significantly improve the 6 interpolated values (RMS values ot 1.369 for back- 
propagation network and 1.029 for radial-basis network). However, with the addition of 10 new points at 
y r =0.025 and 0 075 at A/?=l,2,3,4 and 5, (25-simulation training data set) both the back-propagation 
network and the radial-basis network can accurately capture the overall behavior of the aerodynamic data 
as shown in Figure 43c. The RMS error now is 0. 141 for back-propagation network and 0.106 tor radial- 
basis network. For AR interpolations, the back-propagation network resulted in lower RMS values when 
compared to the RMS values of radial-basis networks (Figure 44). For the 9-point simulation training data, 
the RMS of radial-basis network (RMS=11.12) is quite high when compared to the RMS of back- 
propagation (RMS= 1 . 172) (Figure 44a). For this case, adding 6 new points at AR = 2 and 4 at y r =0, 0.05 and 
0. 1 significantly improves the RMS value for radial-basis (RMS=0.87) as shown in Figure 44b. With the 
addition of 10 new points to 15-simulation data at y ( =0.025 and 0.075 at A/?=l, 2,3,4 and 5, the RMS error 
decreases further to 0.7 for radial-basis networks, and 0.026 for back-propagation (Figure 44c). The results 
indicate that the back-propagation network is quite accurate for small to modest number of data tor the 
cases investigated and it is also more consistent than that of the radial-basis network. However, as indicated 
in Table 32 and Table 33. in terms of computing time or epochs, back- propagation network scales 
unfavorably with respect to the number of data used. In other words, the back-propagation network is 
competitive for modest data size while the radial-basis network is more effective for larger data size. More 
information will be presented when the 2-D airfoil case that involves substantially larger data size is 
discussed. 
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4.6.3 Comparison of Radial-Basis Neural Network and Polynomial-Based Techniques 

For the 3-D wing model, the outputs of the solverb radial-basis NN, along with the results ot the 
polynomial -based technique, are compared for different size of the data. It must be noted that the network 
parameters used to obtain radial-basis network results are 5C= 1.1 75 and error goal=10' 2 . Figure 45 
illustrates the comparison between the NN and polynomial -based outputs based on the 9-simulation 
training data set. For this case, both methods reproduced the original 9-simulation accurately but both tailed 
to predict accurately the interpolation points at y c =0.025 and 0.075 with RMS errors at the test data ot 1 .04 
for both the NN and polynomial-based methods. Furthermore, it is seen that the error estimate of 1. 1 16 of 
Table 3 1 is a gross underestimate. Note that by the time there are 25 data points. Table 31 predicts an error 
of 0.659. The reason for this problem is that RMS error estimates are not reliable when the number of 
coefficients is close to the number of points (7 vs. 9 for this case). In addition, these estimates assume 
random noise and that underlying function is quadratic. Figure 45b shows that adding 6 new r points AR = 2 
and 4 at v ( - = 0, 0.05 and 0. 1 does not help noticeably to improve the 6 interpolated values (RMS values ot 
1.029 for both). However, with the addition of 10 new points aty f =0.025 and 0.0075 at A/?=l,2,3,4 and 5, 
(25-simulation training data set) both the NN and polynomial-based techniques accurately capture the 
overall behavior of the aerodynamic data as shown in Figure 45c. The generalization of the NN with 25- 
simulation is further assessed by comparing additional interpolated values at different y c and AR at 
vc=0.0125 and 0.0875 at M/?=I,2,3,4 and 5. The RMS errors now are 0. 142 for the polynomial and 0.221 
for the NN, which are more in the line with the prediction in Table 31. 

These comparisons illustrate that both neural network and conventional polynomial fitting 
methods do a good job as the number of points is increased. 

5. CONCLUSION AND FUTURE DIRECTIONS 

Recent experiences in utilizing the global optimization methodology, based on polynomial and 
neural network techniques, for aerodynamics and rocket propulsion components are summarized. Global 
optimization methods can utilize the information collected from various sources and by different tools. 
These methods offer multi-criterion optimization, handle the existence of multiple design points and trade- 
offs via insight into the entire design space, can easily perform tasks in parallel, and are often etfective in 
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filtering the noise intrinsic to numerical and experimental data. Another advantage is that these methods do 
not need to calculate the sensitivity of each design variable locally. The global optimization method can be 
particularly effective with either a polynomial-based response surface or a neural network when 
information from different computational, experimental, and analytical sources needs to be assembled. In 
this article, we present recent experiences in utilizing the global optimization methodology for tasks related 
to the preliminary design of a supersonic turbine, multi-criterion design of three different types of injector 
element (shear co-axial, impingement, and swirl co-axial), performance of a low' Reynolds number wing, 
and shape optimization of a turbulent flow diffuser. A successful optimal design often needs to address the 
issues related to the selection of appropriate training data for constructing the global model, employment of 
the statistical and testing tools to identify appropriate global models, existence of multiple design selections 
and related trade-offs, and consideration of noises intrinsic to numerical and experimental data. These 
issues are discussed. It is seen that the global optimization method can naturally take the confidence level 
of the data into account, offers a number of designs with comparable performance, and allows designers to 
make a more informed decision. We have reviewed direct evidences that demonstrate that appropriate 
selection of design points can significantly reduce the number of data required for constructing the global 
model. In particular, while the FCCD approach can be effective with modest number of design variables, 
OA with D-optimal selection criterion seems to be effective when the number of design variables becomes 
higher. Regarding the relative merits between polynomials and neural networks, based on the results 
reviewed, we can make the following summary. 

1. Higher order polynomials usually perform better than lower order polynomials as they have more 
flexibility. However, exceptions have been noticed which demands that appropriate statistical measures be 
taken to determine the best terms to include in an expression. 

2. Both NN and polynomial based RSM can perform comparably for modest data sizes. 

3. Among all the NN configurations, RBNN designed with solve rb seems to be more consistent in 
performance. 

4. Radial basis networks, even when designed efficiently with solverb , tend to have many more 

neurons than a comparable back-propagation with tan-sigmoid or log-sigmoid neurons in the hidden layer. 

The basic reason for this is the fact that the sigmoid neurons can have outputs over a large region of the 
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input space, while radial basis neurons only respond to relatively small regions of the input space. Thus, 
larger input spaces require more radial basis neurons for training. 

5. Configuring a radial basis network often takes less time than that for a back-propagation network 
because the training process for the former is linear in nature. 

6. While the transfer function employed by any neural network is nonlinear in general, RBNN, with 
the combined feature of flexibility and linear regression is more accurate than BPNN, which requires 
solution of nonlinear systems. 

7. The comparisons demonstrate that for this case there are no significant differences between the 
NN and polynomial based RSM. The results of polynomial-based methods, though, suggest that when the 
error is mostly due to modeling rather than noise, the error estimates of the polynomial -based technique can 
be substantially off. 

8. The NN technique has shown the potential of fitting the data much better than the polynomial- 
based technique. However, this was achieved by using the test data to select the parameters like spread 
constant of the NN which appear to greatly affect the predictive accuracy. That is, it was not possible to use 
only the training data to select the best set of parameters. This indicates that because the NNs do not 
provide the statistical information given by polynomial -based methods, using both test data and training 
data is very important in designing the network. 

9. With the large number of points, and the high order polynomial, the statistical predictions of the 
polynomial -based results matched very well the error at the test data. 

10. The neural networks, when trained appropriately, can be used to generate additional data to help 
enhance the data set for constructing the polynomial. Such a combined approach has been demonstrated in 
Shyy et al. [57] for injector design. 

11. The criteria for selecting the data base exhibit significant impact on the efficiency and 
effectiveness of the construction of the response surface. For example, effectiveness of using OA to select 
the database is demonstrated by Papila et al. [50]. 

12. A multi-level approach can be applied to identify the optimal design points with substantially 
higher accuracy. 
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There are a number of outstanding issues that need to be addressed. In the following, we list 
several topics that we consider important for future research. 

1) Is it possible to develop a comprehensive technique by combining NN and polynomial-based RS 
techniques to help reduce the required data size for optimization? 

Specifically, the work done by Rai and Madavan [51], [52] & [53], Madavan et al. [31], and Shyy 
et al. [57] suggests that NN can be effectively used to supplement the existing training data to help 
generating a more accurate polynomial. RBNN may lack satisfactory filtering properties in some cases 
(Papila et al. [49], Vaidyanathan et al. [71]). However, once trained, RBNN can generate additional design 
data to feed the polynomial -based RSM. Polynomials possess the intrinsic filtering capability. The 
evaluation of the nature of the fluctuations from the data generated by RBNN, and the investigation into 
whether polynomials can use the data effectively, is planned. These features have been addressed in this 
article. 

2) What are the kevs to develop a more robust and flexible NN configuration? 

This has been a topic of research for a long time. In this article, an attempt has been made to study 
the training parameters of the different networks used. There are other important issues, which needs to be 
addressed in the future. For example, the possibility of using a more versatile RBNN in terms of a variable 
design parameter, unlike the present situation where the variable has the same value all around the domain, 
needs to be addressed. Ways to determine NN*s performance via statistical tools, especially tor RBNN 
since it employs a linear model to determine the weight associated with each neuron will be investigated. 

3) What is the scaling rule between the number of neurons, and computing time . versus number o} 
input/output variables and the size of the design data? 

There are several rules of thumb for BPNN in the literature (e.g., Greenman [ 14], Carpenter and 
Barthelemy [4], and Fujita [13]) but to our knowledge, no information is provided for RBNN. 
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4) How can one address the need for generating training and testing data most economically and 


effectively? 

The effect of the selection of the design points on accuracy, scaling and performance ot 
polynomial-based RSM has been addressed. The same has yet to be done for NN. 
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7. APPENDIX 
TRAINING DATA 


Table Al. Training Data Sets for 3-D wing model 


Training Data Set #1 
(9-Simulation) 

Training Data Set # 2 
(15-Simulation) 

Training Data Set # 3 
(25-Simulation) 

AR 


WEwam 

AR 

yc 

C l j, 7C d 

AR 

y c 


1 



1 

o 

2.0011 

1 

0.0 


1 

0.05 

4.1224 

1 

0.05 

4.1224 

1 

0.025 

4 

1 


3.6865 

i 


3.6866 

1 

0.05 

4.1224 

3 

0 

5.6398 

2 


4.03 

1 

0.075 

3.99 

3 

.05 

9.6873 

2 

.05 

7.12 

l 

0.1 

3.6866 

3 

0.1 

8.6806 

2 

0.1 

6.34 

2 

0.0 

4.03 

5 

0 

7.9413 

3 

0 

5.6398 

2 

0.025 

7.07 

5 



3 


9.6873 

2 


7.12 

5 

0.1 

12.8951 

3 


8.6806 

2 

0.075 

6.89 




4 

0 

6.92 

2 

0.1 

6.34" 




4 

.05 

11.99 

3 

0.0 

5.6398 




4 

0.1 

10.87 

3 

0.025 

9.64 




5 

0 

7.9414 

3 

0.05 

9.6873 




5 

0.05 

14.0942 

3 

0.075 

9.39 




5 

0.1 

12.8951 

3 

0.1 








4 

0.0 

6.92 







4 

0.025 

11.86 







4 

0.05 

11.99 







4 

0.075 

11.66 







4 

0.1 

10.87 







5 

0.0 

7.9414 







5 

0.025 

13.83 







5 

0.05 

14.0942 







5 

0.075 

13.73 
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Table A2. Test Data Sets for 3-D wing model based on A/? and y c 


Test Set#l for y. 

Test Set#2 for y c 

Test Set#3 for y c 

TestSet#! for AR 

Test Set#2 for AR 

AR 

y c 

AR 

y c 

AR 

y c 

AR 

yc 

AR 

y c 

1 

0.025 

i 

0.025 

1 

0.0125 

2 

0 

2.5 

0 

1 

0.075 

1 

0.075 

1 

0.0875 

2 

0.025 

2.5 

0.025 

3 

0.025 

— 

0.025 

2 

0.0125 

2 

0.05 

2.5 

0.05 

3 

0.075 

2 

0.075 

2 

0.0875 

2 

0.075 

2.5 

0.075 

5 

0.025 

3 

0.025 

3 

0.0125 

2 

0.1 

2.5 

0.1 

5 

0.075 

3 

0.075 

3 

0.0875 

4 

0 

4.5 

0 



4 

0.025 

4 

0.0125 

4 

0.025 

4.5 

0.025 



4 

0.075 

4 

0.0875 

4 

0.05 

4.5 

0.05 



5 

0.025 

5 

0.0125 

4 

0.075 

4.5 

0.075 



5 

0.075 

5 

0.0875 

4 

0.1 

4.5 

0.1 


Table A3. Performance and heat flux responses for O/F = 4 for the shear co-axial injector element. (Table 
A3-5 together contain 45 data points used as the training set) 


O/F 

v/v 0 

Fcomb* til. 

ERE, % 

Q , Btu/in^-sec 

4.0 

4.0 

4.0 

92.9 

0.753 

4.0 

4.0 

5.0 

96.0 

0.753 

4.0 

4.0 

6.0 

97.6 

0.753 

4.0 

4.0 

7.0 

98.6 


4.0 

4.0 

8.0 

99.0 


4.0 

6.0 

4.0 

95.0 

0.928 

4.0 

6.0 

5.0 

97.1 

0.928 

4.0 

6.0 

6.0 

98.5 

0.928 

4.0 

6.0 

7.0 

99.2 

0.928 

4.0 

6.0 

8.0 

99.4 

0.928 

4.0 

8.0 

4.0 

96.6 


4.0 

8.0 

5.0 

98.2 


4.0 

8.0 

6.0 

99.1 

1.10 

4.0 

8.0 

7.0 

99.4 


4.0 

8.0 

8.0 

99.6 
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Table A4. Performance and heat flux responses for O/F = 6 for the shear co-axia! injector element 


O/F 

V/V„ 

Fcombf 

ERE , % 

Q , Btu/in’-sec 

6.0 

4.0 

4.0 

92.9 

0.691 

6.0 

4.0 

5.0 

96.0 

0.691 

6.0 

4.0 

6.0 

97.6 

0.691 

6.0 

4.0 

7.0 

98.6 

0.691 

6.0 

4.0 

8.0 

99.0 

0.691 

6.0 

6.0 

4.0 

95.0 

0.642 1 

6.0 

6.0 

5.0 

97.1 

0.642 

6.0 

6.0 

6.0 

98.5 

0.642 

6.0 

6.0 

7.0 

99.2 

0.642 

6.0 

6.0 

8.0 

99.4 

0.642 

6.0 

8.0 

4.0 

96.6 

0.741 

6.0 

8.0 

5.0 

98.2 

0.741 

6.0 

8.0 

6.0 

99.1 

0.741 

6.0 

8.0 

7.0 

99.4 

0.741 

6.0 

8.0 

8.0 

99.6 

0.741 


Table A5. Performance and heat tlux responses for O/F - 8 for the shear co-axial injector element. 


O/F 

v/v a 

Fcomb' 

ERE, % 

Q , Btu/in 2 -sec 

8.0 

4.0 


92.9 

0.588 

8.0 

4.0 

5.0 

96.0 


8.0 

4,0 

6.0 

97.6 

0.588 

8.0 

4.0 


98.6 

0.588 

8.0 

4.0 

8.0 

99.0 

0.588 

8.0 

6.0 

4,0 

95.0 


8.0 

6.0 


97.1 


8.0 

6.0 


98.5 


8.0 

6.0 

7 0 

99.2 

0.512 

8.0 

6.0 

8.0 

99.4 

0.512 

8.0 

8.0 


96.6 


8.0 

8.0 

5.0 

98.2 

0.493 

8.0 

8.0 

6.0 

99.1 

0.493 

8.0 

8.0 

7.0 

99.4 

0.493 

8.0 

8.0 

8.0 

99.6 

0.493 
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Table A6. Data used to test the polynomials and NN for the shear co-axial injector element. (The table 
contains 20 data points used as the testing set) 


OIF 

V/V„ 

E com b) til* 

ERE , % 

Q , Btu/in 2 -sec 

4.0 

5.0 

4.0 

94.4 

0.812 

4.0 

5.0 

5.0 

96.9 

0.812 

4.0 

5.0 

6.0 

98.1 

0.812 

4.0 

5.0 

7.0 

99.1 


4.0 

5.0 

8.0 

99.4 


4.0 

7.0 

4.0 

96.0 

1.014 

4.0 

7.0 

5.0 

97.9 

1.014 

4.0 

7.0 

6.0 

98.8 

1.014 

4.0 

7.0 

7.0 

99.4 

1.014 

4.0 

7.0 

8.0 

99.6 

1.014 

6.0 

5.0 

4.0 

94.4 

0.642 

6.0 

5.0 

5.0 

96.9 

0.642 

6.0 

5.0 

6.0 

98.1 

0.642 

6.0 

5.0 

■hbh 

99.1 


6.0 

5.0 


99.4 

0.642 

6.0 

7.0 

4.0 

96.0 

0.691 

6.0 



97.9 

0.691 

6.0 

7.0 

6.0 

98.8 

0.691 

6.0 

7.0 

7.0 

99.4 

0.691 

6.0 

7.0 

8.0 

99.6 

0.691 


Table A7. Propellant momentum ratio as a function of propellant pressure drops: shear co-axial 


injector element. 



APo 

A P f 

200 

180 

160 

150 

140 

120 

100 

200 

1.49 

1.42 

1.33 

1.30 

1.25 

1.16 

1.06 

*— 1 
oo 
o 

1.57 

— r^o 

1.41 

1.37 

1.32 

1.22 

1.11 

160 

1.67 

1.59 

1.50 

1.45 

1.40 

1.30 


150 

1.73 

1.64 

1.54 

1.49 

1.44 

1.34 


140 

1.79 

1.70 

1.60 

1.55 

1.50 

1.39 



1.93 

1.83 

1.72 

1.67 

1.61 

1.50 


100 

2.11 

2.00 

1.89 

1.83 

1.77 

1.64 

1.49 
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Table A8. Design data for a shear co-axial injector element with AP t? and AP. = ZOOpsi. 


AP 0 

AP, 

mm 

a 

ERE 

Q- 

H unpin Ee 

w rtI 

C„, 

200 

200 


15 

NA 

0.85 

0.84 

0.923 

1.083 

200 

200 

2 

20 

85 

0.85 

0.62 

0.923 


200 

200 

2 

30 

92.8 

0.85 

IHsEeHI 



200 

200 

2 

45 

95.4 

BSEEH 


| 

it#! 

200 

200 

2 

50 

95.8 

HEEEH 

0.19 


MEM 










200 

200 

4 

15 

91 

0.85 


1 


200 

200 

4 

20 

95.2 



1 

1.083 

200 

200 

4 

30 


0.85 

0.39 

1 

■EHM 

200 

200 

4 

45 


0.85 

0.23 

1 

■sum 

200 

200 

4 

50 

98.4 

0.85 

0.19 

1 

1 .083 I 
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Figure 2. Schematic of a simple neuron model 


Input Layer 


Hidden Layers 


Output Layer 







75 





Figure 6. Face centered composite designs (FCCD) tor 3 design variables; x h .vj, and xj 



Figure 7. The two phases of the optimization process, where Phase 1 deals with data processing/generating 

and Phase 2 deals with optimization. 
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Figure 8. Desirability function for various weight factors, s. 



Figure 9. Schematic of shear co-axial injector element 
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Figure 10. Schematic of impingement injector element 
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Figure 16. Comparison of a for different NN designed with solverb for shear co-axial gas-gas injector (a) 

ERE (Cr) and (b) Q ( Btu/in : -sec). 
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(a) (t>) 

Figure 18. Predictive capabilities of various models tor ERE of the shear co-axial injector, (a) Polynomial 

(b) NN. 



(a) (b) 

Figure 19. Predictive capabilities of various models for Q of the shear co-axial injector, (a) Polynomial (b) 
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Figure 20. Comparison between cubic and quadratic response surface for case 3 of (a) ERE. (b) Q for the 
shear co-axial injector. (NN and injector model data are the same points in the graph. Quadratic and cubic 


are predicted by polynomial -based RSM) 



Figure 21. Assessment of predictive capabilities of (a) polynomial based method, (b) RBNN for ERE of 
the shear co-axial injector. (15 training/mapping points 10 test points) 
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Figure 22. Assessment of performance of cubic and quadratic response surfaces tor case 1 of (a) ERE ot 
the shear co-axial injector. (25 training/mapping points) (b) Q. ( 13 training/mapping points) { Enhanced set 
includes Injector model data and RBNN predicted data. Quadratic and cubic are predicted by polynomial- 


based RSM} 


Comparison txivwsm cubic and quaAulic response surfaces 
Case 1 . s * t . t = 10. Lcomh = 7 
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Figure 23. Assessment of performance of cubic and quadratic response surfaces for case I of (a) ERE for 
the shear co-axial injector. (25 training/mapping points) (b) £>.(13 training/mapping points) { Enhanced set 
includes Injector model data and RBNN predicted data. Quadratic and cubic are predicted by polynomial- 


based RSM} 
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Figure 24. Performance of various models for ERE of the shear co-axial injector, fa) RS’s, (b) NN’s. 



(a) (b) 

Figure 25. Performance of various models for Q of the shear co-axial injector, (a) RS’s, (b) NN’s. 
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Figure 26. Performance and weight trends for swirl and impinging elements. 



Figure 27. Heat flux and cost trends for swirl and impinging elements. 






Figure 28. Effect of the number of turbine stage on optimum design parameters; D, RPXt, and A ann and 
optimum output parameters; rj, VV’ and Apay calculated for Apay-based optimization (All geometric design 



Figure 29. Comparison of crfor different NN designed with solverbe for single-stage supersonic turbine (a) 


rj i r/ c) and (b) \V \ lbs). 
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rj (%) and (b) W (lbs). 
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Figure 31. Comparison of error goal vs. number of epochs and erfor networks trained with solverb for 
single-stage supersonic turbine (a) rj (7c) and (b) W (lbs). 
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Effect of Constraints on Optimum RPM 
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Figure 32. Effect due to presence (case 1) or lack of constraints (case 2) on design variables, (a) optimum 
diameter, D (in ), (b) optimum RPM, (c) optimum annulus Area, A ann (in. 2 ), (d) optimum vane axial chord, 
Cv (in.), (e) optimum blade axial chord, C b (in.), (0 optimum stage reaction, K r (%). 
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(1: No constraints: 2: Constraints 
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(e) 

Figure 33. Effect due to presence (case 1) or lack of constraints (case 2) on objective functions, (a) 
Optimum Efficiency, rj, (b) Optimum Weight. W (lbs), (c) Optimum pitch speed, V plIch (in. /sec), (d) 
Optimum (A, ln „ x RPM), AN 2 (in : *rpm 2 ). (e) Optimum Incremental Payload. Apay (lbs). 
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Nounalzed Value of DV n (-1.+1) 



(a) Original Design Space (b) Refined Design Space 


Figure 34. Comparison of the design variables for optimization based on payload increment (Apay) using 
1990-data (FCCD) and 249-data (OA) for both original design space and refined design Space ( DV#l : D, 
DV#2: RPM , DV#3: A ann , DV#4: h h DY#5: c,/, DV#6: c, :> DV#7: c (ti , DV#8: c h2 , DV#9: sr It DVttlO: sr 2l 
and DVUl I: w /} ). Both designs are satisfactory, demonstrating that there exist multiple optimum designs. 



Figure 35. Effect of spread constant \sc) on training rms-error (cr a ) of (a) rf and (b) W tor preliminary 
design of 2-stage turbine for original design space using 1990-training data tor solverbe RBNN 
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(a) (W 


Figure 36. Effect of spread constant ( 5 c) on training rms-error (a a ) of (a) ij and (b) W for preliminary 


design of 2-stage turbine for original design space using 249-training data for solverbe RBNN 



(a) for r| (Apay optimization) (b) lor W (Apay optimization) 

Figure 37. Effect of spread constant (sc) on training rms-error (a a ) of for preliminary design of 2-stage 
turbine for refined design space using 249-training data (OA) for solverbe, RBNN. 
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249 training data (OA), and 78 testing data in both ori 

are normalized by t 


9 - 




Figure 40. Normalized revalues along a straight line in the design space. The results are for two different 
differencing schemes and two different residual levels used as convergence criterion. (UDS: upwind 
differencing scheme and SUDS: second-order upwind differencing scheme) 
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Figure 41. Comparison of optimum wall shapes using polynomial and B-spline representations, 
respectively. 





i f 

:/ 


i 


Figure 42. Contour plot of response surface for diffuser design. Solid circle indicates the optimal region. 

The hatched part of the feasible space comprises designs with performance within \% of the optimal. 
Corresponding shapes are indicated to the right The results indicate that multiple design points meet the 


2oal. 
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(a) 9-Simulation Training (b) 15-Simulation Training (c) 25-Simulation Training 


Figure 43. Comparison of radial basis network with back-propagation network results for 3-D wing model 
(for w interpolation) (Design parameters: sc=l 175, error goal=0. 1 for radial-basis, and #of neurons=4. 


error goal =0.01 for back-propagation) 



(a) 9-Simulation Training <b) 15-Simulation Training (c) 25-Simulation Training 


Figure 44. Comparison of radial basis network with back-propagation network results for 3-D wing model 
(for AR interpolation) (Design parameters: yc= 1.175, error goal=0.1 tor radial-basis, and #ot neurons— 4, 


error goal =0.0 1 for back-propagation) 



(a) 9-Simulation Training (b> 15-Simulation Training (c) 25-Simulation Training 


Figure 45. Comparison of radial basis NN results with polynomials for 3-D wing model (Design 
parameters: sc=\A15 . error goal=0. 1 for radial-basis networks) 
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10. TABLES 


Table 1. Comparison of NN and polynomial- based response surface (RS) techniques 



NN- based 
RS\1 

Polynomial- 
based RSM 

Comments 

Computational Effort 
and Cost 

Disadvantage 

Advantage 

Finding the weights associated with the neurons is a 
non-linear regression process tor all of the NN 
types other than RBNN. Whereas, finding the 
polynomial coefficients requires solution of a linear 
set of equations. 

The cost increases if the regression process is non- 
linear which makes NN’s other than RBNN more 
expensive than polynomials. 

Noise 

Disadvantage 

Advantage 

Ability of filtering noise from experimental data is 
possible with polynomial-based RSM. However, it 
the number of neurons used to design the NN is not 
the same as the data, then, by definition, filtering is 
also possible for NN-based RSM. 

Handling complex 

functions 

Advantage 

Disadvantage 

NN’s are more suitable for multi-dimensional 
interpolation of data that lack structure since they 
are much more flexible in functional form 
especially when dealing with design in the context 
of unsteady flows, partial and/or complete data sets. 


Table 2. Literature review on NN and polynomial-based RS techniques comparison 


Authors 

No of 
Data 

No of 
Input 

No of 
Output 

NN Type 
(2- layer) 

Activation 

Function 

No of 
Neurons 

Polynomial 

Degree 

Carpenter & 

Barthelemy 

[4] 

36 

2 

1 

BPNN 

Sigmoid 

1, 2, 4 

l 5 ' to 4 tn 

961 

2 

1 

3, 5. 7 

2 nd to 5 ,n 1 

81 

4 

1 

1,2, 3 

1“ to 2 nd 1 

300 

15 

1 

2,4, 6, 8, 10 

I s1 to 2 nd 

Madavan et al. 

- 


nH 

BPNN 

Sigmoid 

15 & 7 


1 st to 2 nd 

[31] 



1 

(3- layer) 





Nikolaidis et 

400 

50 

i 

BPNN 

Sigmoid 

a of NN is insensitive 

^nd 

al. [39] 






to no of neurons 


Papila et al. 

9 

2 

i 

RBNN 

radbas 

8, 9 

4 


[49] 

15 

2 

i 

& 

& 

12, 15 

4 

gSBEHMI 


25 

2 

i 

BPNN 

Sigmoid 

20, 25 

4 




2 

i 



253, 255 

- 




2 

i 



765 

- 


Rai and Madavan 

3 & 5 

1 

i 

BPNN 

Sigmoid 

1 & 2 


[52] 

27 

3 

i 

(3- layer) 


7 & 3 

■HSESIH 


. 

15 

i 



- 


Shyy et al. [57] 

45 

3 

2 

RBNN 

radbas 

42 and 45 


Vaidyanathan et 

45 

3 

2 

RBNN 

radbas 

42 and 45 



al. [71] 

76 

6 

2 





11 


98 

























Table 3. Design variables and design space tor single-, 2- and 3-Stage turbines 


(All geometric design variables are normalized by the baseline values) 


SINGLE-STAGE 

2-STAGE 

3-STAGE 


Variable 

Lower 

Limit 

Upper 

Limit 

Variable 

Lower 

Limit 

Upper 

Limit 

Variable 

Lower 

Limit 

Upper 

Limit 

Mean 

Diameter, D 

0.50 

1.50 

Mean 

Diameter, D 

0.50 

1.50 

Mean 

Diameter, D 

0.50 

1.50 

Speed, RPM 

0.70 




1.30 

Speed, RPM 

0.70 

1.30 

Blade 

Annulus 

Area, 

0.70 

1.30 


0.70 

1.30 

Blade 
Annulus 
Area, A*™ 

0.70 

1.30 

Vane Axial 
Chord, c v 

0.39 

1.71 

1 st Blade 
Height (% of 
Exit Blade), 
hi 

0.90 

1.50 ; 

1 st Blade 
Height (% 
of Exit 
Blade) , h { 

0.90 

1.50 

Blade Axial 
Chord, c b 

0.26 

1.14 

1 st Vane 
AxtaJ Chord, 

Cvl 

0.39 

1.71 

1 st Vane 
Axial 
Chord, Cvi 

0.39 

1.71 

Stage 

Reaction, sr 

0.0% 

50% 

1 st Blade 
Axial Chord, 
CbL 

0.26 

1. 14 

I s ' Blade 
Axial 
Chord, Cbi 

0.26 

1.14 



2 nd Vane 
Axial Chord, 

C V 1 

0.21 

1.41 

2 ,ld Vane 
Axial 
Chord, c V 2 

0.21 

1.41 


2 nd Blade 
Axial Chord, 

Cb2 

0.17 

1.13 

2 nd Blade 
Axial 
Chord, c b2 

0.17 

1.13 

I st Stage 
Reaction, sr f 

0.0% 

50% 

3 rd Vane 
Axial 
Chord, c v3 

0.21 

1.41 

2 nd Stage 
Reaction, sr> 

0.0% 

50% 

3 rd Blade 
Axial 
Chord, c b3 

0.17 

1.13 

1 st Work 
Fraction, wf L 

50% 

85% 

I st Stage 
Reaction, sri 

0.0% 

50% 



2 nd Stage 
Reaction, sr 2 

0.0%) 

50% 


3 rd Stage 
Reaction, sr 3 

0.0% 

50% 

1 st Work 
Fraction, wf[ 

40% 

80% 

2 nd Work 
Fraction, wf 2 

30% 

10% 
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Table 4. Different cubic polynomials for ERE . (Dependent variables: V/V a and L rombl 15 training points, 10 


test points). 


Model 

# 

Coefficient = 0 

Terms 

Removed 

Terms Included 

<7a (%) 

cr(%) 

1 

V/Vo 3 



0.09 

0.21 

2 

V/Vo 3 

V/Vo~L com h 


0.08 

0.21 

3 


V/v7l~. vyvo’ 


0.08 

0.21 

4 


V/V7L~. V/Vo' 

£-comb 

0.09 

0.21 

5 


V/v7L~. V/Vo’ 

lZ7. V/V,;L, ,~- 

0.09 

0.21 

6 


V/v7l~, V/Vo 1 ~ 

E C omb < V/Va E C omb > 
V/VoL,„ mh - 

0.10 

0.21 


Table 5, Different cubic polynomials for Q. (Dependent variables: O/F and V/V 0 , 9 training points, 4 test 

P° intS )- * 


Model 

# 

Coefficient = 0 

r / 

Terms 

Removed 

Terms Included 

<Ja (%) 

<T(%) 

1 

vyvv/, o/f j 



5.58 

2.23 

2 

O/F 1 

VA'o 2 


5.58 

2.09 

3 


V/Vo 7 , O/F 3 


5.58 

2.09 

4 


V/W. O/F 7 


5.58 

2.23 

5 


V/Vo*. O/F'. V/Var 


3.96 

2.09 

6 


V/Vo'. Q/F\ V/Vo' 

V/Vo- O/F 

5.58 

2.09 


Table 6. Neural Network architectures used to design the model for shear co-axial injector element, {sc - 
spread constant } 


Scheme 

# of 

Layers 

t — : ’■ 

# of neurons in j # of neurons in 
the hidden j the output layer 

layer > 

Error goal aimed for during training 



ERE 

Q ERE 

Q 

ERE 

Q 

RBNN 
(Solve rbe) 

2 

15 

9 : i 

1 

0.0 { sc = 3.2 5) 

0.0 {jc= 1.20} 

RBNN {Solverb) 

2 

14 

8 1 

l 

0.001 {«' = 1 05 1 

0.00 1 {.sc =1.05} 

BPNN 

2 

8 

4 i 1 

l 

0.01 

0.01 


Table 7. RMS error in predicting the values of the objective function by various schemes for the shear co- 


axial injector element. 


Scheme 

a for ERE (%) 

a for Q ( % ) 

RBNN ( Solverbe ) 

0.20 

1 40 

RBNN ( Solverb ) 

0.13 

r 1.53 

BPNN 

0.18 

0.83 

Partial Cubic RS 

0.21 

2.23 ' 

Quadratic RS 

0.28 

3.49 
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Table 8. Optimum values obtained with cubic and quadratic for case 1. (Constraints. 4 < O/F < 6, 4 < 
Vf /y o < 6, and U < 7) (Values in the parenthesis are the exact resp onse of the injector model] 


Vff 

v 0 - u, a 

tu Mromb 

Cubic 

' 

~ r - 


Quadratic 



Were 

(s) 

w q 

(0 

O/F 

v/v 0 

Leo mb 

ERE 

Q 

O/F 

Vf/V„ 

Lcomb 

ERE 

Q 

1 

10 

6.0 

5.41 

7.0 

99.02 

(99.00) 

0.664 

(0.654) 

6.0 

6.00 

7.0 

99.17 

(99.20) 

0.669 

(0.642) 

l 

1 

6.0 

6.00 

7.0 

99.15 

0.669 

6.0 

6.00 

7.0 

99.17 

0.669 

10 

1 

6.0 

6.00 

7.0 

99.15 

0.669 

6.0 

6.00 

7.0 

99.17 

0.669 


Table 9. Optimum values obtained with cubic and quadratic for case 2. (Constraints. 4 < O/F < 6 t 5 < 
V/V 0 < 7, and Lcomh ^ 7) (Values in the parenthesis are the exact response of t he injector model} — 



v 0 b /, ai 

— 

1U Lcomb —'M valuer in .. 

Cubic 

Quadratic 










Q 

O/F 

v/v. 

i 

o 

J 

- 

ERE 

!■ 






99.02 

(99.00) 

0.664 

(0.654) 

6.0 

6.52 

7.0 

99.31 

(99.10) 

U.6S4 

(0.716) 

■ 

m 




99.21 

(99.20) 

0.674 

(0.691) 

6.0 


7.0 

99.42 

(99.30) 

0.702 

(0.728) 






99.32 


6.0 


7.0 

99.42 



Table 10. Optimum values obtained with cubic and quadratic for case 3. (Constraints: 4 < O/F < 6, 6 < 

V r /Vo< 8, and Lc omb <7) 



— 

v v o — «•> *-comD - 

Cubic 

Quadratic 




Were 

(s) 

Wq 

(t) 

O/F 

v/v. 

Lcomb 

ERE 

Q 

O/F 

V/V. 

Lcomb 

ERE 

Q 

l 

10 

6.0 

6.00 

7.0 

99.15 

0.669 

6.0 

6.52 

7.0 

99.31 

0.684 

l 

1 

6.0 

6.34 

7.0 

99 21 

0.674 

6.0 

8.00 

7.0 

99.67 

0.753 

10 

1 

6.0 

8.00 

7.0 

99.42 

0.728 

6.0 

8.00 

7.0 

99.67 

0.753 
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Table 11. The RMS error in the prediction of ERE and Q for different values ot spread constant. The error 


goal used for Solve rb is 0.001. 



Solverbe 

.W. 

Solverbe 

Solverb 

Solverb 

Solverb 

sc 

RMS error 
(ERE) 

RMS error 

(Q) 

RMS error 
(ERE) 

RMS error 

(0) 

No of neurons 

0.50 

1.493 

0.179 

1.733 

0.287 

44 

0.75 

0.745 

0.135 

0.675 

0.135 

44 

1.00 

0.152 

0.022 

0.153 

0.017 

42 

1.05 

0.190 

0.011 

0.128 

0.012 

44 

r i.25 

0.316 

0.010 

0.267 

0.022 

44 

1.50 

0.336 

0.022 

0.309 

0.030 

44 

1.75 

0.369 

0.022 

0.310 

0.021 

44 

2.00 

0.308 

0.016 

0.296 

0.019 

41 

2.25 

0.279 

0.020 

1.846 

0.045 

43 

2.50 

0.325 

0.017 

0.744 

0.025 

43 


Table 12. Optimum values obtained with cubic and quadratic for case l (enhanced data set). (Constraints. 




H- KJ! I — P v O — '•'I M:omD — ' / 

Cubic 

Quadratic 



Were 

(s) 

Wq 

(t) 

O/F 

v/v 0 

Lcomb 

ERE 

Q 

O/F 

WVo 

Leo mb 

ERE 

Q 

l 

10 

6.0 

5.54 

7.0 

99.02 

(98.90) 

0.654 

(0.658) 

6.0 

5.01 

7.0 

98.96 

(98.70) 

0.644 

(0.664) 

1 

1 

6.0 

6.00 

7.0 

99.12 

0.658 

6.0 

6.00 

7.0 

99.25 

0.658 

10 

1 

6.0 

6.00 

7.0 

99.12 

0.658 

6.0 

6.00 

7.0 

99.25 

0.658 


Table 13. Optimum values obtained with cubic and quadratic for case 2 (enhanced data set ). (Constraints. 




Cubic 

Quadratic 



WfcKE 

(S) 

W Q 

(t) 

O/F 

V,/V 0 

Lcomb 

ERE 

Q 

O/F 

V./V,, 

Lcomb 

ERE 

Q 

1 

10 

6.0 

5.54 

7.0 

99.02 

(98.90) 

0.654 

(0.658) 

6.0 

5.01 

7.0 

98.96 

(98.70) 

0.644 

(0.664) 

1 

1 

6.0 

6.33 

7.0 

99.18 

(99.10) 

0.663 

(0.666) 

6.0 

6.04 

7.0 

99.26 

(99.20) 

0.659 
(0 642) 

10 

j 

1 

6.0 

7.00 

7.0 

99.30 

0.681 

6.0 

7.00 

7.0 

99.46 

0.693 
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Table 14. Optimum values obtained with cubic and quadratic for case 3 (enhanced data set). (Constraints. 

4 < Q/F < 6, 6 < V/Vq < 8, and Uomb < 7 ) {Compare with Table 10J 




Cubic 

Quadratic 




W ER e 

(s) 

Wq 

(t) 

O/F 

V/Vo 

Lcomb 

ERE 

Q 

O/F 

V/V 0 

L«:omb 

ERE 

Q 

1 

10 

6.0 

6.00 

7.0 

99.12 

0.658 

6.0 

6.00 

7.0 

99.25^ 

0.658 

1 

1 

6.0 

6.33 

7.0 

99.19 

0.663 

6.0 

6.04 

7.0 

99 26 

0.659 

10 

1 

6.0 

8.00 

7.0 

99.42 

0.725 

6.0 

7.95 

7.0 

99.57 

0.746 


Table 15. Optimal Solutions for fixed values of V/V 0 and given range ot O/F and L aim b obtained with NN 
and RSM schemes for the shear co-axial injector element. (Constraints: 4 < O/F < 8, 4 < L tomb < 7) (error 

^ C7~\ 


v/v 0 

given in 

Scheme 

pareninesib iui 

O/F 

CaclJ pltuiCLivm to 

in. 

ERE, % 

O , Btu/in 2 -sec 

4 

RBNN (Solve rbe) 

8.0 

7.0 

98.60 (0.00) 

0 588 (0.00) 

RBNN ( Solverb ) 

8.0 

7.0 

98.60 (0.00) 

0 588 (0.00) 

BPNN 

8.0 

6.9 

98.64(0.14) 

0.578 (1.70) 

Partial Cubic RS 

8.0 

7.0 

98.61 (0.01) 

0.595 (1.19) 1 

Quadratic RS 

8.0 

7.0 

98.67 (0.07) 

0.591 (0.51) 

Model 

8.0 

7.0 

98.60 

0.588 

Model 

8.0 

6.9 

98.50 

0.588 

6 

RBNN (Solve rbe) 

8.0 

7.0 

99.20 (0.00) 

0.512 (0.00) 

RBNN ( Solverb ) 

8.0 

7.0 

99.20 (0.00) 

0.512(0.00) 

BPNN 

8.0 

7.0 

99.18 (0.02) 

0.513 (0.20) 

Partial Cubic RS 

8.0 

7.0 

99.15 (0.05) 

0.499 (2.54) 

Quadratic RS 

8.0 

7.0 

99.17 (0 03) 

0.531 (3.71) 

Model 

8.0 

7.0 

99.20 

0.512 

8 

RBNN ( Solverbe ) 

8.0 

7.0 

99.40 (0.00) 

0.493 (0.00) 

RBNN ( Solverb ) 

8.0 

7.0 

99.40 (0.00) 

0.493 (0.00) 

BPNN 

S.O 

7.0 

99.41 (0.01) 

0 500(1.42) 

Partial Cubic RS 

8.0 

7.0 

99.42 (0.02) 

0 500(1.42) 

Quadratic RS 

8.0 

7.0 

99.67 (0.27) 

0 471 (4.46) 

Model 

S.O 

7.0 

99.40 i 0.493 
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Table 16. Effect of each variable on the optimization of impingement co-axial injector element -optimal 
designs for original constraints & equal weights 


Independent 

Variable 

Constraints 

Results Case 1 

Results Case 2 

E — 

Results Case 3 

Results Case 4 

AP n 

100-200 

183 

183 

179 

100 

u q 

APr 

100-200 

100 1 

132 

149 

100 

I 

2-8 

8.0 

8.0 ^ 

6.6 ^ 

6.5 

a 

15-50 

33.1 

' 18.9 

22.3 

24.0 

Dependent 

Variable 

Desirability 

Limits 

ERE & Q» 

ERE, Q w > 

H impinge 

ERE, Q w , 

Himpinge? ^^rel 

ERE, Qw^ 

Himpioge* 

Crel 

ERE 

95.0-99.9 

99.9 

98.3 

98.0 

98.0 

Ow 

0.7-1 .3 

0.74 

0.76 

H 0.79 

0.86 

HuBDinff 

0.2-1 .0 

— 

0.75 

0.61 

0.63 

w„, 

0.9-1 .2 

— 

— 

1.1 

1.1 

Crtl 

0.7-1. 1 

— 

— 


0.93 


Table 17. Effect of emphasizing & life & performance issues on the optimization of impingement co-axial 
injector element —optimal designs for original constraints and moditied weights 


in 

Independent 

Variable 

jeciur cicmciu - 

Constraints 

Results 

Casel 

Constraints 

Results 
Case 2 

Constraints 

Results 
Case 3 

AP 0 

100-200 

100 

100-200 

158 

100-200 

100 

APr 

100-200 

100 1 

100-200 

100 

100-200 

137 


2-8 

6.5 

2-8 

7.7 

2-8 

5.2 

ct 

15-50 

24.0 

15-50 

15.0 

15-50 

36.0 

Dependent 

Variable 

Baseline 

Variable 

Weight 


Life 

Variable 

Weight 


Thrust fW eight 
Variable 
Weight 


ERE 

1 

98.0 

1 

96.7 

5 

99.1 

Ow 

1 

0.86 

5 

0.75 

r i 

0.95 

V* 

H- 

1 

0.63 

5 

0.94 

i 

0.32 

iMnioinee 

VV reI 

1 

1.10 

1 

1.14 

5 

1 .05 

C rt , 

1 

0.93 

1 

0.97 

1 

0.95 
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Table 18. Effects of realizing the last increments of performance & weight on the optimization of 


impingement 

Independent 

Variable 

LU-tiAiai 

Original 

Constraints 

Results 

Casel 

- 

Modified 

AP 

Constraints 

Results 
Case 2 

Modified 

AP 

Constraints 

Results 
Case 3 

AP 0 

100-200 

100 

150-200 

150 

150-200 

150 

AP f 

100-200 

137 

150-200 

150 

150-200 

200 


2-8 

5.2 

2-8 

5.4 

2-8 

4.4 

a 

15-50 

36.0 

15-50 

35.6 

15-50 

44.8 

Dependent 

Variable 

Variable 

Weight 

(5:1) 


Variable 

Weight 

(5:1) 


Variable 

Weight 

(100:1) 


ERE 

5 

99.1 

5 

99.0 

10 

99.1 

Qw 

1 

0.95 

1 

0.84 

0.1 

0.95 

Ejnujlnft 

1 

0.32 

1 

0.31 

0.1 

0.21 

w rel 

5 

1.05 

5 

1.05 

10 

1.01 

c„, 

1 

0.95 

1 

1.00 

0.1 

1.07 


Table 19. Effect of each variable on the optimization of swirl co-axial injector element —optimal designs 

for original constraints and equal weights. T 


Independent 

Variable 

Constraints 

Results 
Case 1 

Results 
Case 2 



Results 
Case 3 

Results 
Case 4 

Results 
Case 5 

AP 0 

100-200 

200 

200 

200 

200 

104 

APf 

20-200 

41 

41 

42 

47 

20 

i 

2-8 

7.2 

7.2 

7.6 

3.2 

3.4 

© 

30-90 

81 

81 

37 

47 

44 ^ 

Dependent 

Variable 

Desirability 

Limits 

ERE 

ERE & Q„ 

ERE, Qw, 
Qinj 

ERE, Q„, 
Qinj* Wrel 

ERE, Q w , 
Qinj* W rel , 
Crel 

ERE 

92.3-99.0 

98.5 

98.5 

97.2 

96.0 

95.7 

Qw 

0.596- 

0.647 

0.596 

0.596 

0.596 

0.596 

0.596 

J&JI 

Qini 

6.95-36.59 

26.8 

26.8 

9.1 

12.0 

10.5 



W ret 

0.900- 

1.154 

1.13 

1.13 

1.14 

0.97 

0.98 

C„i 

0.73-1 .42 

0.98 

0.98 

0.81 

0.84 

0.76 
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Table 20. Effect of emphasizing life and performance issues on the optimization of swirl co-axial injector 


element. 


Independent 

Variable 

Constraints 

Results 

Baseline 

Constraints 

Results 
Case 1 

Constraints 

Results 
Case 2 

AP 0 

100-200 

104 

100-200 

200 

100-200 

200 



AP f 

20-200 

20 

20-200 

32 

20-200 

44 


2-8 

3.4 

2-8 

3.6 

2-8 

2.9 

0 

30-90 

44.0 

30-90 

30.0 

30-90 

72.0 

Dependent 

Variable 

Baseline 

Variable 

Weight 


Life 

Variable 

Weight 


Thrust/Weight 

Variable 

Weight 


ERE 

1 

95.7 

1 

95.3 

10 

96.7 

Ow 

1 

0.596 

5 

0.596 

1 

0.596 



Ow 

1 

10.5 

10 

6.9 

1 

22.6 



w rel 

1 

0.98 

1 

0.99 

2 

0.96 

C„, 

1 

0.76 

1 

0.79 

1 

0.94 


Table 21. Optimization summary for l, 2 and 3-stage turbine with response surface in original design 



r r 

Original Design Space 

Refined Design Space 

non, Wop, Apay„p, 

F|ooi W op t Apay 0P t 

1 -stage 

Apay 2-stage 

3-stage 

0.77 0.73 -0.21 

1.10 1.05 0.11 

1.24 1.62 0.14 

0.77 0.73 -0.21 

1.13 1.04 0.15 

1.20 1 54 0.11 


Table 22. Values of <7 ( , and crfor different response surfaces of rj and W for the supersonic turbine. 


Scheme 

(7 a for T](%) 

crfor r\ (%) 

<x for W (%) 

<rfor W(%) 

Quadratic RS 

2.51 

0.90 

0.82 

1.27 

Reduced Cubic RS 

1.95 

1.03 

0.40 

1.22 


Table 23. Neural Network architectures used to design the models for rf and W of the supersonic turbine. 


Scheme 

# of 

Lavers 

# of neurons in the 
hidden layer 

# of neurons in the 
output layer 

Error goal ain 
training 

led for during 



n 

W 

i 

W 

n 

VV 

RBNN (Solve rbe) 

2 

76 

76 

l 

1 

0.0 

{.sc = 9.50} 

0.0 

{$c = 9.45 } 

RBNN (Solve rb) 

2 

37 

! 

75 

1 

1 

0.001 

{.sc = 6.50} 

0.001 

{sc = 8.35 } 

BPNN 

2 

5 

60 

1 

l 

0.001 

0.001 
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Table 24. RMS error in predicting the values of rj and Why various schemes for the supersonic turbine. 


Scheme 

a for r}(%) 

0 for W (%) 

RBNN ( Solverbe ) 
RBNN ( Solverb ) 

1.25 

0.29 

1.10 

1.10 

BPNN 

Reduced Cubic RS 

0.78 

1.03 

2.56 

1.22 


Table 25 Optimal solutions with constraints on V pMh and AN 2 for a supersonic turbine, (error given in 
2 each peed,™™ is in %) (V pLIC i, = 1. 148 and AN 2 = 1 . 132 in all .he cases) (All , he variables 

are normalized by their respective baseline values) 


Scheme 


RBNN 
1 (Solve rbe) 


D , in. 


0.972 


0.972 


RPM 


1.181 


1.181 


0.811 


0.811 


C„ in. 


1.443 


1.443 


C by in. 


0.836 


0.836 


Kn 


0.0 


0.0 


0.810 

(5.80) 


0.766 


W , lbs 


0.636 

(0.74) 


0.641 


Apay 


lbs 


-0.139 

(29.80) 


- 0.197 


RBNN 
(Solve rb) 


Meanline 


BPNN 


Meanline 


Reduced 
Cubic RS 
Meanline 


0.999 


0.999 


1.024 


1.024 


0.903 


0.903 


1.149 


1.149 

1.121 


0.857 

0.857 


1.483 


0.792 


0.0 


0.785 

(1.75) 


0.653 

(0.17) 


1.483 


0.792 


0.0 


0.772 


0.654 


0.901 


1.168 


1.143 


1.121 


0.901 


1.168 


1.143 


0.0 

oo” 


0.793 

(2.49) 


0.608 

(8.63) 


0.772 


0.666 


1.272 


1.272 


0.700 


0.700 


1.706 


1.706 


0.871 


0.0 


0.758 

(1.50) 


0.871 


0.0 


0.746 


0.591 

( 2 . 10 ) 


0.604 


-0.177 
(9 16) 


- 0.194 


-0.153 

(21.49) 


- 0.195 


-0.194 

(8.40) 


- 0.211 


Table 26. The quality of the second-order response surface obtained for 77 , W and Apay of 2-stage turbine 
for 1990-data (face centered criterion) and 249-data (orthogonal arrays) (Mean values ot rj, VVand J P QV are 




ri 

W 

Apay 


R 7 

0.995 

0.996 

0.995 


Ra 2 

0.994 

0.996 

0.995 

1990-data 

Rms- error 

1.31% 

2.56% 

9.58% 


Mean 

0.78 

0.86 

-0.24 


R" 

0.995 

0.998 

0.994 


Ra 2 

0.992 

0.997 

0.992 

249-data 

Rms- error 

2.128% 

0.826%) 

20 . 68 % 


Mean 

0.89 

0.92 

- 0. 1 1 
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Table 27. Testing of the second-order response surface obtained for r/and W of 2-stage turbine for 1990 

° v 1 /AA \ IQ Ant** 


data 

# of design points 

llllCll'Jii; anu if/ 

# of test data 

a for rj (%) 

a for W(%) 

249 

78 

P 1.65 

0.96 

1990 

78 

1.67 

1.21 


Table 28. Testing 

the RBNN and second-order polynomial response surface obtained tor rj and W for 
nreliminarv design of 2-stage turbine (original design space) - 

Number of 
training 
data 

Number of sc 

test data 

a for r\ (%) 
using RBNN 
(Solverbe) 

a for n (%) using a for W(%) 

Polynomial -based using RBNN 

RSM (Solverbe) 

o for W(%) using 
Polynomial -based 
RSM 

249 

1990 

78 

78 

4.3 

3.2 

1.365 

2.263 

1.648 1-305 

1.672 1-557 

0.959 

1.214 


Table 29. Backward elimination procedure for polynomial -based RSM in five variables. 


Terms 

Min |to| 

No. |t 0 | <2.0 

- ^ 

R,r 

Comments 

21 

0.05 1 

15 1 

' 0.922 "1 

0.811 

^ A 3 

?0 

0.23 

14 

0.922 

0.823 

Removed 

™ 

19 

0.45 

12 

0.922 

0.834 

Removed y 4 

18 

0.53 

9 

0.921 

0.841 

Removed V/V4 

17 

0.97 

8 

0.919 

0.848 

Removed v>v< 
^ T -1 7 

*16 

1.22 

h 6 

0.915 

0.848 

Removed y 2 

15 

1.57 

5 

0.909 

0.844 

Removed y/yj 


Table 30. Comparison between CFD-solutions and polynomial -based RSM-predictions. 


y, 

y 2 

Y\ 

Yi 

Ys 

F 

F (full) 

F (reduced) 

0 5 

0.5 

0.5 

0.5 

0.5 

0.7171 

0.7148 

0.7126 

1.0 

p; 0.5 

0.0 

0,5 

1.0 

0.7174 

0.7210 

0.7174 

0.25 

0.75 1 

0.25 

0.75 

0.25 

0.7148 

0.7185 

0.7162 

0.0 

0.5 

1.0 

0.5 

0.0 

0.6943 

0.7333 

0.7283 
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Table 31. Predicted RMS error, a for different polynomial models for 3-D wing model: 9-simulation, 
15-simulation, and 25-simulation data sets. (The shaded models indicate the best fit) 


Model 

No 

MODEL 

a for 
9 data 

a for 
15 data 

a for 
25 data 

1 

Ci A R‘+C 2 A R+CjA Ry c +c 4 y c +c 5 y<f+c 6 


■illM 


2 

C|AR“+C2AR+C3ARyc + c -<yc + c syc~+ c 6 + C 7 AR 3 


KMi 


3 

CiAR 2 +C 2 AR+c 3 ARy c +c 4 yc +c 5 yc“+c 6 + c?yc J 


HSjjgWg 

0.5524 ! 

4 

CiAR^+CjAR+CiARyc +c 4 yc +Csyc~+ C 6 + C7ARyc“ 


0.0738 


5 

CiAR 2 +c 2 AR+c 3 ARy c +c 4 y c +c 5 y c ‘+c 6 +c 7 yc J +c 8 ARy c " 

- 

■■KM 

0.3207 

6 

CiAR"+C 2 AR+C}ARyc +c 4 yc +Csyc 2 +C 6 

+c 7 yc 3 +cgARy c 2 +C9ycAR 2 

- 

■ 


7 

ciAR 2 +C 2 AR+c 3 ARy c +Ciy c +c 5 yc 2 +C 6 
+ c 7 y r 3 +cgARyc'+ C 9 AR 

- 

- 

0.6961 

8 

CiAR z +C 2 AR+c 3 ARy c +c 4 y c +c 5 y c '+c 6 
+c 7 y c 3 +cgARyc 2 +C9ycAR 2 + C 10 AR 3 

- 

- 

0.3350 

9 

C]AR+C 2 ARyc +C3yc +C4yc 2 +Cs +C6yc 3 +C7ARyc 2 

- 

- 

0.4248 

10 

CiAR+c 2 ARy c +c 3 y c +c 4 y c 2 +Cs -K: 6 AR 2 yc 

- 

- 

0.8044 

11 

c j AR”+C 2 A R+c 3 A Rye +c 4 yc +Csyc 2 +Ct 
+c 7 y c 3 +c g ARy c 2 +C 9 yc 4 

- 

- 

0.2383 

12 

c 1 AR / +C 2 AR+c 3 ARy c +c 4 y c +c 5 yc 2 +C 6 
+c 7 y c 3 +CgARyc 2 +C9yc 4 +CioARyc 3 

' 


0.1073 


Table 32. Training history of radial -basis networks with Solve rb for 3-D wing model 


NN 

No. 

#of 

Simulations 

# of Neurons 

# of Epochs 

Steady State 
Error 

Spread 

Constant 

Error Goal 

1 

9 

8 

7 

10 16 

1.175 

10 2 

2 

15 

12 

11 

10* 4 

1.175 

10 2 

3 

25 

20 

19 

io 3- 

1.175 

10 2 


Table 33. Training history of back-propagation networks with Trainlm for 3-D wing model 


NN 

No. 

#of 

Simulations 

# of Neurons 

# of Epochs 

Steady State 
Error 

Error Goal 

1 

9 

4 

23 

4.5x10 4 

10* 2 

2 

15 

4 

12 

8.5xlO ' 3 

10' 2 

3 

25 

4 

105 

9.96x10 3 

10 " 
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